/*
**   maxlik12.e  -  A Tobit model with Bayesian inference
**
*/

library maxlik,pgraph;
#include pgraph.ext;
#include maxlik.ext;
graphset;
maxset;

proc lpr(x,z);
    local s,m,u;
    s = x[4];
    if s <= 0;
       retp(error(0));
    endif;
    m = z[.,2:4]*x[1:3,.];
    u = z[.,1] ./= 0;
    retp(u.*lnpdfn2(z[.,1]-m,s) + (1-u).*(ln(cdfnc(m/sqrt(s)))) );
endp;


rndseed 13731;

start = { 1, 1, 1, 1 };
__title = "tobit example";
output file = maxlik12.out reset;

{x0,f0,g0,cov0,ret0} = maxlik("tobit",0,&lpr,start);
call MAXPrt(x0,f0,g0,cov0,ret0);

cl0 = MAXTlimits(x0,cov0);
call MAXClprt(x0,f0,g0,cl0,ret0);

_max_NumSample = 200;
_max_BootFname = "mbayes";
__output = 10;


proc prior(b);   /*  unit normal prior  */
   retp(prodc(pdfn(b)));
endp;
_max_PriorProc = &prior;

__title = "Weighted Bayesian Bootstrap";
{x1,f1,g1,cov1,ret1} = MAXBayes("tobit",0,&lpr,x0);

call MAXPrt(x1,f1,g1,cov1,ret1);

cl1 = MAXBlimits(_max_BootFname);
call MAXClprt(x1,f1,g1,cl1,ret1);

printdos "\n  MAXBayes has created the GAUSS data file `mbayes'\n";
printdos "  containing the weighted bootstrap coefficients.  You may\n";
printdos "  want to delete this file after the program has been run\n\n";
"press any key for kernel density plots...";
call keyw;


call MAXDensity(_max_BootFname,0);

output off;