samer@0: % mscaler - Dynamic additive and multiplicative normalisation in Matlab
samer@0: %
samer@0: % mscaler ::
samer@0: %    model(real) ~'model for observations',
samer@0: %    options {
samer@0: %       scale   :: nonneg /1  ~'initial scale factor';
samer@0: %       offset  :: real   /0  ~'initial offset';
samer@0: %       scale_rate  :: nonneg/0.02 ~'scale adaptation rate';
samer@0: %       offset_rate :: nonneg/1e-7 ~'offset adaptation rate';
samer@0: %       batch   :: natural/4  ~'batch size for updates'
samer@0: %    }
samer@0: % -> arrow({[[N]]},{[[N]]},mscaler_state).
samer@0: 
samer@0: function o=mscaler(model,varargin)
samer@37: 	opts=options('scale',1,'offset',0, ...
samer@0: 		'scale_rate',0.02,'offset_rate',1e-7,'batch',4, ...
samer@0: 		'nargout', 1, ...
samer@0: 		varargin{:});
samer@0: 
samer@0: 	rates=[opts.scale_rate;opts.offset_rate];
samer@0: 	batch=opts.batch;
samer@0: 	score=scorefn(model);
samer@0: 	z3 = [0;0;0];
samer@0: 
samer@0: 	if (opts.batch==1)
samer@0: 		o=loop(@update,@(s)[opts.scale;opts.offset]);
samer@0: 	else
samer@0: 		o=loop(@update_batched,@(s){z3,[opts.scale;opts.offset]});
samer@0: 	end
samer@0: 
samer@0: 	function ss=stats(y,phi), 
samer@0: 		n=size(y,1); 
samer@0: 		ss=[n;sum(y.*phi)-n;sum(phi)]; 
samer@0: 	end
samer@0: 
samer@0: 	function [y,phi]=infer(params,x)
samer@0: 		y = (x-params(2))/params(1);
samer@0: 		phi = score(y);
samer@0: 	end
samer@0: 
samer@0: 	function params=updparams(params,stat1,stat2)
samer@0: 		params = [	params(1)*exp(rates(1)*stat1); 
samer@0: 						params(2)+rates(2)*params(1)*stat2 ];
samer@0: 	end
samer@0: 
samer@0: 	function [y,state]=update(x,state)
samer@0: 		[y,phi] = infer(state,x);
samer@0: 		state = updparams(state, mean(y(:).*phi(:))-1, mean(phi(:)));
samer@0: 	end
samer@0: 
samer@0: 	function [y,state]=update_batched(x,state)
samer@0: 		params=state{2};
samer@0: 		[y,phi] = infer(params,x);
samer@0: 		ss = state{1} + stats(y(:),phi(:));
samer@0: 		if (ss(1)<batch) state{1}=ss;
samer@0: 		else state = {z3, updparams(params,ss(2)/ss(1),ss(3)/ss(1))};
samer@0: 		end
samer@0: 	end
samer@0: end