samer@0: % afitgauss - Fit Gaussian to stream of input vectors
samer@0: %
samer@0: % afitgauss :: 
samer@0: %    nonneg~'adaptation rate',
samer@0: %    options {
samer@0: %       zeromean :: bool/0 ~'assume input is zero mean?'
samer@0: %    }
samer@0: % -> arrow( {[[N]]}, {struct { mu::[[N]]; cov:[[N,N]] }}, stats_struct(N)).
samer@0: %
samer@0: % stats_struct(N) ::= struct { 
samer@0: %    ag  :: nonneg; 
samer@0: %    amu :: [[N]]; 
samer@0: %    aco ::[[N,N]];
samer@0: % }.
samer@0: 
samer@0: function o=afitgauss(L,varargin)
samer@0: 
samer@37: 	opts=options('zeromean',0,varargin{:});
samer@0: 
samer@0: 	GMM=gmm;
samer@0: 	%K=L/(L+1);
samer@0: 	kappa=L*(L-1); % precision of diffusion process for effective memory of L
samer@0: 
samer@0: 	M0=[];
samer@0: 	o=loop(ifx(opts.zeromean,@update_z,@update),@init);
samer@0: 
samer@0: 	function s0=init(size_in)
samer@0: 		N=size_in(1);
samer@0: 		M0=GMM.ss_init(GMM.blank(1,N));
samer@0: 		s0=foldl( @(s,x)GMM.ss_accum(s,[x,-x],[0.5,0.5]), M0, mat2cell(eye(N),N,ones(1,N)));
samer@0: 	end
samer@0: 
samer@0: 	function [y,s]=update(x,s)
samer@0: 		% factor to reduce weight of evidence so far
samer@0: 		% corresponds roughly to *addition* of variance
samer@0: 		% in the exponential family diffusion scheme of things
samer@0: 		K=kappa/(kappa+s.ag); 
samer@0: 		s.ag  = K*s.ag + size(x,2);
samer@0: 		s.amu = K*s.amu + sum(x,2);
samer@0: 		s.aco = K*s.aco + x*x';
samer@0: 		y=mstep(s,M0);
samer@0: 	end
samer@0: 
samer@0: 	function [y,s]=update_z(x,s)
samer@0: 		% factor to reduce weight of evidence so far
samer@0: 		K=kappa/(kappa+s.ag); 
samer@0: 		s.ag  = K*s.ag + size(x,2);
samer@0: 		s.aco = K*s.aco + x*x';
samer@0: 		y=mstep_z(s,M0);
samer@0: 	end
samer@0: 
samer@0: 	function M=mstep(S,M)
samer@0: 		M.mu=S.amu/S.ag;
samer@0: 		M.cov=msym(S.aco/S.ag - M.mu*M.mu'); 
samer@0: 	end
samer@0: 
samer@0: 	function M=mstep_z(S,M)
samer@0: 		M.mu=zeros(size(S.aco,1),1);
samer@0: 		M.cov=S.aco/S.ag; 
samer@0: 	end
samer@0: 
samer@0: 	function A=mysm(X), A=(X+X')/2; end
samer@0: end