proc calprobs1(h,startt,maxt);
/* turns a hazard rate into a distribution */
/* same as calprobs */
local probs, alive, i,hh;

probs=zeros(maxt,2);
hh=h[.,1]+h[.,2];
hh=1-hh[startt+1:startt+maxt];
alive=cumprodc(hh);
alive=1|alive[1:maxt-1];

probs[.,1]=alive.*h[1+startt:maxt+startt,1];
probs[.,2]=alive.*h[1+startt:maxt+startt,2];

retp(probs);
endp;

proc calprobs(h,startt,maxt);
/* turns a hazard rate into a distribution */
local probs, alive, i;

probs=zeros(maxt,2);
alive=1;
for i (1,maxt,1);
probs[i,1]=alive*h[i+startt,1];
probs[i,2]=alive*h[i+startt,2];
alive=alive*(1-h[i+startt,1]-h[i+startt,2]);
endfor;

retp(probs);
endp;

proc(3)=petersen(h1,startt,maxt);
local prue,probs1,probs2,pp1,pp2,ff,fmax,fhat,tt,ptrue,i,fmin, smin,smax,it;

probs1=calprobs1(h1,startt,maxt);
pp1=probs1[.,1]~probs1[.,2];


fmin=zeros(maxt*100,1); fmax=fmin+1;
tt=seqa(0.00051,0.01,maxt*100);
for i(1,maxt*100,1);
it=int(tt[i]);
if (it.ge 1).and (it.le maxt); fmin[i]=sumc(pp1[1:it,2]); endif;
if (it.lt maxt); fmax[i] =sumc(pp1[1:it+1,2])+sumc(pp1[1:it+1,1]); endif;
endfor;
retp(tt,1-fmax,1-fmin);
endp;

proc(1)=plopeter(tt,smi,sma,tit1,tit2);
local h1, a1,a2;

xy(tt,sma~smi); 
retp(0);
endp;






library pgraph;
graphinit;


begwind;
window(2,2,0);

tt=seqa(45,1,36);

setwind(1);
_pltype=6|4|6|4; 
_plwidth=2;
_pcolor=1;
_pnumht=.25;
_ptitlht=0.25;
fonts("simplex complex microb simgrma");
/* setup legend */
_plegctl = { 1 3 47 0.051 };
_plegstr = " 1970 lower bound\0 2000 lower bound\0 1970 upper bound\0 2000 upper bound";



scale(45~80,0~1); 
xtics(45,80,10,0);
ytics(0,1 ,0.1,0);
_plctrl=0;

load h1=haz_70wm;
load h2=haz_80wm;
load h3=haz_90wm;
load h4=haz_00wm;

startt=45; maxt=35;
{tt,s70l,s70h}=petersen(h1,startt,maxt);
{tt,s00l,s00h}=petersen(h4,startt,maxt);
tit1="bounds";
tit2="White Males";
plopeter(tt+45,s70l~s00l,s70h~s00h,tit1,tit2);


endwind;