Code.m

Image processing

xmean=[120e3 160e3];

covx=[0.3.3];

ygroup=[10;01];

xr=eye(length(xmean));

save Praxmean xr covx ygroup

%%

for i =1 :size(ygroup,1)

yc=-MBeta*ygroup(i,:);

[ysamples]=Buildysamples_normal(yc,1);

parfor j = 1:size(ysamples,1)

[Fs(i,j),Gy(i,j)]=CallFLAC_RSM(j,ysamples(j,:)]

end

[RSCCoef(i,:)]=BuildRS(Gy(i,:),ysample]

[yd(i,:),beta(i,:),flag(i,:)]=FORM(RScoef(i,:),yc]

end

%%

ym=zeros(1,length(xmean));

ysd=ones(length(xmean),1);

yr=xr;

yc1=yd(1,:);

yc2=yd(2,:);

Cy=(ysd*transpose(ysd)).*yr;

N=100;

pf1=normcdf(-beta(1),0,1);

pf2=normcdf(-beta(2),0,1);

a1=pf1/(pf1+pf2);

a2=1-a1;

parfor i=1:N

u=rand(1,1);

if u<w1

yc=yc1;

else

yc=yc2;

end

y=mvnrnd(yc,Cy,1);

w(i,1)=(mvnpdf(y,ym,Cy))/(a1*mvnpdf(y,yc1,Cy)+a2*mvnpdf(y,yc2,Cy));

[Fsm(1,i),Gsm(1,i)]=CallFLAC_IS(i,y)

if Fsm(1,i)>1

1(i,1)=0;

else

1(i,1)=1;

end

end

pf_m=mean(1.*w)

pf_std=std(1.*w)/sqrt(N)

cov= pf_std/pf_m

sysbeta=-norminv(pf_m,0,1)

save RSTcovx=[0.3.3];ygroup=[10;01];xr=eye(length(xmean));save Praxmean xr covx ygroup%%for i =1 :size(ygroup,1) yc=-MBeta*ygroup(i,:); [ysamples]=Buildysamples_normal(yc,1); parfor j = 1:size(ysamples,1) [Fs(i,j),Gy(i,j)]=CallFLAC_RSM(j,ysamples(j,:)] end [RSCCoef(i,:)]=BuildRS(Gy(i,:),ysample] [yd(i,:),beta(i,:),flag(i,:)]=FORM(RScoef(i,:),yc]end%%ym=zeros(1,length(xmean));ysd=ones(length(xmean),1);yr=xr;yc1=yd(1,:);yc2=yd(2,:);Cy=(ysd*transpose(ysd)).*yr;N=100;pf1=normcdf(-beta(1),0,1);pf2=normcdf(-beta(2),0,1);a1=pf1/(pf1+pf2);a2=1-a1;parfor i=1:N u=rand(1,1); if u<w1 yc=yc1; else yc=yc2; end y=mvnrnd(yc,Cy,1); w(i,1)=(mvnpdf(y,ym,Cy))/(a1*mvnpdf(y,yc1,Cy)+a2*mvnpdf(y,yc2,Cy)); [Fsm(1,i),Gsm(1,i)]=CallFLAC_IS(i,y) if Fsm(1,i)>1 1(i,1)=0; else 1(i,1)=1; endendpf_m=mean(1.*w)pf_std=std(1.*w)/sqrt(N)cov= pf_std/pf_msysbeta=-norminv(pf_m,0,1)save RST