%enter the following
cin = 'in';      %control file
out = 'out';     %output log file


%do not write anything below
LD=50000;
LDPAR=300;
LVAL=100;
s=' ';

%open the control file and log file
fin=fopen(cin);
if fin==-1; disp('failure to open file %s'); disp(cin); quit; end;    
fout=fopen(out,'w');
if fout==-1; disp('failure to open file'); disp(out); quit; end;

%read the control file
fprintf(fout,'%10sPrior density ratio class robustness\n\n',s);
nd=fscanf(fin,'%d',1); 
nburn=fscanf(fin,'%d',1);
ifunc=fscanf(fin,'%d',1);   
if (nd<1) | (nd>LVAL);
   fprintf(fout,'**Number of d values %4d not in range 1 to %4d\n',kern,LVAL); 
   quit;
end;
dpar=plistr(nd,fin,fout);
fprintf(fout,'  Number of burn-in iterations:%10.d\n',nburn);
fprintf(fout,'          Function of interest:%10.d\n',ifunc); 
hfile=fscanf(fin,'%s',1);

%open the input file and output file
in=fopen(hfile);
if in==-1; fprintf(fout, 'failure to open file %s', hfile); quit; end;



%open simulation file and read
niter=fscanf(in,'%d',1);   errcheck(niter,fout,0);
npars=fscanf(in,'%d',1);    errcheck(npars,fout,0);

%error check
if (nburn<0) | (nburn>niter);
   fprintf(fout,'**Burn-in iterations %4d not in range 0 to niter',nburn,niter);
   quit;
end;
if (ifunc<1) | (ifunc>npars);
   fprintf(fout,'**Function of interest %4d not in range 1 to npars',ifunc,npars);
   quit;
end;

nuse=niter-nburn;
nskip=floor((nuse-1)/LD)+1;
for iburn=1:nburn;
   [jter, wtlog, pri, pdat, par]=readrec(npars,in,fout);
end;

%pass through file to get weights and values
m=0;
for iter=1:nskip:nuse;
   m=m+1;
   [jter, wtlog, pri, pdat, par]=readrec(npars,in,fout);
   if iter==1; wtmax=wtlog; end;
   if (wtlog>wtmax) wtmax=wtlog; end;
   wt(m)=wtlog;
   g(m)=par(ifunc);
end; %for iter

%convert from log weights to weigths and find moments
mtot=m;
for m=1:mtot;
   wt(m)=exp(wt(m)-wtmax);
end; %for m

[amu,asd,qup]=rob1(mtot,wt,g,nd,dpar);
g=g*(-1);
[aa,asd,qdown]=rob1(mtot,wt,g,nd,dpar);
qdown=qdown*(-1);

global ak;
for id=1:nd;
   ak=exp(dpar(id));
   f0=0;
   f=fsolve('robeq',f0);
   beta(id)=(ak-1)*phi(f)/(ak+(1-ak)*cumnord(f));
   fprintf(fout,'%15.5e\n',beta(id));
end; %for id
fprintf(fout,'%15sRobustness analysis of function%3d\n\n',s,ifunc);
fprintf(fout,' Posterior mean%15.6e   Stan dev%15.6e\n\n',amu,asd);
fprintf(fout,'%24sExact',s);
fprintf(fout,'            DeRobertis & Hartigan asymptotic\n');
fprintf(fout,'%9sd    Lower bound    Upper bound',s);
fprintf(fout,'    Lower bound    Upper bound\n');
for i=1:nd
   fprintf(fout,'%10.6f%15.6e%15.6e%15.6e%15.6e\n',dpar(i),qdown(i),qup(i), amu-asd*beta(i),amu+asd*beta(i));    
end;
fclose(in);
fclose(fout);
quit;