0001 function [beta,g,P2,P0]=gmem(x,f,beta0,varargin);
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 TOL=.005;
0015 options=strvcat(varargin{:});
0016 sf=sum(f);
0017 f=f(:)/sf;
0018 x=x(:);
0019 A=beta0(:,1);
0020 m=beta0(:,2);
0021 r=beta0(:,3);
0022 Ni=length(A);
0023 Nx=length(x);
0024 dx=(x(end)-x(1))/(Nx-1);
0025 minr=dx.^2;
0026
0027 beta=beta0; err=inf;
0028 while err>TOL,
0029
0030
0031 xm=(x(:,ones(1,Ni)).'-m(:,ones(1,Nx))).^2;
0032 P0=exp(-xm./(2*r(:,ones(1,Nx))))./sqrt(2*pi*r(:,ones(1,Nx)))*dx;
0033 P1=A(:,ones(1,Nx)).*P0;
0034
0035
0036 g=sum(P1,1);
0037 P2=P1./(eps+g(ones(1,Ni),:));
0038
0039
0040
0041 A=eps+P2*f;
0042
0043 m=(P2*(f.*x))./A;
0044
0045 if strmatch('comvar',options),
0046 r=sum((P2.*xm)*f)*ones(Ni,1);
0047 else
0048 r=((P2.*xm)*f)./A;
0049 end
0050 r=max(minr,r);
0051
0052 betanew=[A,m,r];
0053 err=abs(betanew-beta)./(eps+beta); err=max(err(:));
0054 beta=betanew;
0055 end
0056
0057 beta(:,1)=beta(:,1)*sf;
0058 g=g*sf;
0059