view arrows/stats/afitgauss.m @ 61:eff6bddf82e3 tip

Finally implemented perceptual brightness thing.

% afitgauss - Fit Gaussian to stream of input vectors
%
% afitgauss :: 
%    nonneg~'adaptation rate',
%    options {
%       zeromean :: bool/0 ~'assume input is zero mean?'
%    }
% -> arrow( {[[N]]}, {struct { mu::[[N]]; cov:[[N,N]] }}, stats_struct(N)).
%
% stats_struct(N) ::= struct { 
%    ag  :: nonneg; 
%    amu :: [[N]]; 
%    aco ::[[N,N]];
% }.

function o=afitgauss(L,varargin)

	opts=options('zeromean',0,varargin{:});

	GMM=gmm;
	%K=L/(L+1);
	kappa=L*(L-1); % precision of diffusion process for effective memory of L

	M0=[];
	o=loop(ifx(opts.zeromean,@update_z,@update),@init);

	function s0=init(size_in)
		N=size_in(1);
		M0=GMM.ss_init(GMM.blank(1,N));
		s0=foldl( @(s,x)GMM.ss_accum(s,[x,-x],[0.5,0.5]), M0, mat2cell(eye(N),N,ones(1,N)));
	end

	function [y,s]=update(x,s)
		% factor to reduce weight of evidence so far
		% corresponds roughly to *addition* of variance
		% in the exponential family diffusion scheme of things
		K=kappa/(kappa+s.ag); 
		s.ag  = K*s.ag + size(x,2);
		s.amu = K*s.amu + sum(x,2);
		s.aco = K*s.aco + x*x';
		y=mstep(s,M0);
	end

	function [y,s]=update_z(x,s)
		% factor to reduce weight of evidence so far
		K=kappa/(kappa+s.ag); 
		s.ag  = K*s.ag + size(x,2);
		s.aco = K*s.aco + x*x';
		y=mstep_z(s,M0);
	end

	function M=mstep(S,M)
		M.mu=S.amu/S.ag;
		M.cov=msym(S.aco/S.ag - M.mu*M.mu'); 
	end

	function M=mstep_z(S,M)
		M.mu=zeros(size(S.aco,1),1);
		M.cov=S.aco/S.ag; 
	end

	function A=mysm(X), A=(X+X')/2; end
end