#include <oxstd.h>
#import "bprobit"

main()
{
    decl probitobj, ct = 100, x, y, crep = 100, irep,
        ires, cfailed, params, parses;

    probitobj = new Bprobit();
    probitobj.Create(1, 1, 1, ct, 1);   // create database
    probitobj.Deterministic(FALSE);     // create constant

    x = ranu(ct, 1);            // fixed during experiment
    y = zeros(ct, 1); // 0 as yet, created in replications
    probitobj.Append(x ~ y, {"x", "y"});

    probitobj.Select(Y_VAR, { "y",0,0 } );    // formulate
    probitobj.Select(X_VAR, { "Constant",0,0, "x",0,0 } );
    probitobj.SetSelSample(-1, 1, -1, 1);   // full sample

    probitobj.SetPrint(FALSE);   // no intermediate output
    params = parses = <>;
    for (irep = cfailed = 0; irep < crep; )
    {
        y = 1 + x + rann(ct, 1);  // create new y variable
        y = y .< 1 .? 0 .: 1;             // make into 0,1
        probitobj.Renew(y, {"y"});  // replace in database

        probitobj.ClearModel();     // force re-estimation
        ires = probitobj.Estimate();

        if (probitobj.GetModelStatus() != MS_ESTIMATED)
        {
            ++cfailed;             // count no of failures
            continue;      // failed: reject and try again
        }
        params ~= probitobj.GetPar();  // store parameters
        parses ~= probitobj.GetStdErr(); // and std errors
        ++irep;                        // next replication
    }

    println("No of successful replications: ",
        crep, " (", cfailed, " failed)");
    println("Sample size: ", ct);
    println("estimated parameters", 
        "%c", {"mean-par", "sd-par", "mean-se", "sd-se"},
        meanr(params) ~ sqrt(varr(params)) ~
        meanr(parses) ~ sqrt(varr(parses)));

    delete probitobj;
}
