#include <oxstd.h>
#import <maximize>

decl g_mY;                                  // global data
decl g_mX;                                  // global data

fProbit(const vP, const adFunc, const avScore,
    const amHessian)
{
    decl prob = probn(g_mX * vP);   // vP is column vector
    decl tail = 1 - prob;

    adFunc[0] = double(
        meanc(g_mY .* log(prob) + (1-g_mY) .* log(tail)));

    if (avScore)                   // if !0: compute score
    {
        decl weight = (g_mY - prob) .* densn(g_mX * vP)
            ./ (prob .* tail);	
        avScore[0] = meanc(weight .* g_mX)';// need column
		print("avScore[0]", meanc(weight .* g_mX)');// it's indeed a column vector, 3x1
    }

return 1;                           // 1 indicates success
}


main()
{
    decl vp, dfunc, ir, mhess, cn;

    print("Probit example 2, run on ", date(), ".\n\n");

    decl mx = loadmat("data/finney.in7");

    g_mY = mx[][0];       // dependent variable: 0,1 dummy
    g_mX = 1 ~ mx[][3:4]; // regressors: 1, Lrate, Lvolume
    cn = rows(g_mY);
    delete mx;

    vp = <-0.465; 0.842; 1.439>;        // starting values
    MaxControl(-1, 1, 1);          // print each iteration
                                               // maximize
    ir = MaxBFGS(fProbit, &vp, &dfunc, 0, FALSE);

    print("\n", MaxConvergenceMsg(ir),
        " using analytical derivatives",
        "\nFunction value = ", dfunc * cn,
        "; parameters:", vp);

                  // if converged: compute standard errors
    if (ir == MAX_CONV || ir == MAX_WEAK_CONV)
    {
        if (Num2Derivative(fProbit, vp, &mhess))
        {
            mhess = -invert(mhess) / cn;
            print("standard errors:", sqrt(diagonal(mhess)));
        }
    }
}