/* enter the following */
nd=3;       /* number of d values */
nburn=100;  /* burn in iterations */
ifunc=3;     /*function of interest */
dpar={0.69 1 2}; /* d values */

hfile="../../../models/uvr1/sim"; /* posterior simulation file */
outfile="out"; /* output log file name */



      /* do not write anything below */
let LDPAR=300; 
let LD=50000;
let LVAL=100;
tabfile="table";

output file=^outfile reset; 

/* check user info */
print "          Prior density ratio class robustness";
print;
if ((nd<1) or (nd > LVAL));
   format /rdn 4,0; 
   print "Number of d values " nd "not in range 1 to" LVAL; end;
endif;
for i (1,nd,1);
   if (dpar[i]<=0);
      format /rdn 4,0; 
      print "Element "  i "of list ";;
      format /len 15,7;
      print "not positive"; end;
   endif;
endfor;

/* write user info to log file */
format /rdn 10,0;
print "  Number of burn-in iterations: " nburn;
print "          Function of interest: " ifunc;    


load in[]=^hfile;
niter=in[1];
npars=in[2];
if not (rows(in) == (niter*(4+npars) + 2));
   print "
          Error in simulation input file line 1 
          or unexpected end of file reached at line " (niter*npars +2) "
          or some of the data is not numeric ";
   end;
endif;
if ((nburn<0) or (nburn > niter));
   format /rdn 4,0; 
   print "Burn-in iterations " nburn "not in range 1 to " niter; end;
endif;
if ((ifunc<1) or (ifunc > npars));
   format /rdn 4,0; 
   print "Function of interest" ifunc "not in range 1 to " npars; end;
endif;

nuse=niter-nburn;
nskip=floor((nuse-1)/LD+1);
in=in[3:rows(in)];
in=reshape(in, niter, npars+4);
jter=in[.,1];
mtot=floor((nuse-1)/nskip+1);
gw=zeros(mtot,2);
/*gw=selif(in,jter.>nburn .and nskip.*floor((jter+1)./nskip).==(jter+1));*/
gw[.,2]=selif(in[.,2],jter.>nburn .and nskip.*floor((jter+1)./nskip).==(jter+1));
gw[.,1]=selif(in[.,(4+ifunc)],jter.>nburn .and nskip.*floor((jter+1)./nskip).==(jter+1));
wtmax=maxc(gw[.,2]);
gw[.,2]=exp(gw[.,2]-wtmax);

q=zeros(1,nd);
amu=0; asd=0;
call rob1(gw);
qup=q;
gw[.,1]=-gw[.,1];
call rob1(gw);
qdown=-q;
amu=-amu;


load gparin[]=^tabfile;
gparin=reshape(gparin, 100, 2);
beta=zeros(1,nd);
for i (1,nd,1);
   if (dpar[i] <= gparin[1,1]);
      f= dpar[i]*gparin[1,2]/gparin[1,1];
      beta[i]=f;
   elseif ((dpar[i]>gparin[1,1]) and (dpar[i] <= gparin[100,1]));
      fu=selif(gparin,gparin[.,1].>=dpar[i]);
      fl=selif(gparin,gparin[.,1].<dpar[i]);
      ll=rows(fl);
      f=fl[ll,2] + (dpar[i]-fl[ll,1])*(fu[1,2]-fl[ll,2])/(fu[1,1]-fl[ll,1]);
      beta[i]=f;
   endif;
endfor;



format /rdn 4,0;
print "               Robustness analysis of function" ifunc;
print;
format /len 15,6;
print " Posterior mean  " amu "Stan dev  " asd;
print;
print "                        Exact";;
print "           DeRobertis & Hartigan asymptotic";
print "         d     Lower bound    Upper bound";;
print "    Lower bound    Upper bound";

for i (1,nd,1);
   format /rd 10,6;
   print dpar[i];;
   format /len 15,6;
   if (dpar[i] <= gparin[100,1]);
      print qdown[i] qup[i] amu-asd*beta[i] amu+asd*beta[i];
   else;
      print qdown[i] qup[i] "d value too large";
   endif;
endfor;
   
end;