function [L,t,p1,w]=DAmixture(x,tnew,pnew) % modification of Data algorithm (see Liu p 134 ) applied to a normal variance mixture with % 2 components % test using e.g. [t,w,p1]=DAmixture(normmixrnd(50,[0 0],[1 3],[2/3 1/3])); % for the actual data augmentation algorithm, use DAmix.m t=tnew+.001; p1=pnew+.001; %tnew= preliminary estimates of the standard deviation of the two components % p1=preliminary estimates of the probability of the first normal % population L=[]; f1=normpdf(x,0,t(1)); f2=normpdf(x,0,t(2)); w=(p1*f1)./(p1*f1+(1-p1)*f2); while (sum(abs(t-tnew))+abs(p1-pnew))>.0000001 t=tnew; p1=pnew; tnew=sqrt([sum(w.*x.^2)/sum(w) sum((1-w).*x.^2)/sum(1-w)]); f1=normpdf(x,0,tnew(1)); f2=normpdf(x,0,tnew(2)); w=(p1*f1)./(p1*f1+(1-p1)*f2); pnew=mean(w); L=[L sum(w.*log(f1)+(1-w).*log(f2))]; % the vector of values of the loglikelhood end