#include <oxstd.h>

#include "bprobit.h"


Bprobit::Bprobit()
{
    Modelbase();                   // initialize base class
}

Bprobit::fProbit(const vP, const adFunc, const avScore, const amHess)
{
    decl prob = probn(m_mX * vP);
    decl tail = 1 - prob;

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

    if (avScore)
    {
        decl weight = (m_mY - prob) .* densn(m_mX * vP)
            ./ (prob .* tail);
        avScore[0] = meanc(weight .* m_mX)';
    }

return 1;                            // 1 indicates success
}

Bprobit::InitPar()
{
	Modelbase::InitPar();   // first call Modelbase version
	SetParCount(m_cX);
	SetPar(zeros(m_cX, 1));              // start from zero
    return TRUE;
}
Bprobit::DoEstimation(vStart)
{
    m_iResult = MaxBFGS(fProbit, &vStart, &m_dLogLik, 0, FALSE);
       // m_iResult has return code, vStart will have new params
    m_dLogLik *= m_cT;             // change to unscaled log-lik
    return {vStart, "BFGS", FALSE};       // three return values
}
Bprobit::Covar()
{
    decl mhess;
    if (Num2Derivative(fProbit, m_vPar, &mhess))
    {
        m_mCovar = -invert(mhess) / m_cT;
    }
}   
