Function __mlogistic_trans_compute

src/modules/regress/clustered_errors.cpp:394–472

void __mlogistic_trans_compute (MutableClusteredState& state,
                               const MappedColumnVector& x, const double& y)
{
	int numCategories = state.numCategories -1;
	ColumnVector yVec(numCategories);
    yVec.fill(0);

    //Pivot around the reference category
    if (y > state.refCategory) {
        yVec((int)y - 1) = 1;
    } else if (y < state.refCategory) {
        yVec((int)y) = 1;
    }

//    if ((int)y != 0) {
 //       yVec(((int)y) - 1) = 1;
 //   }

	/*
    Compute the parameter vector (the 'pi' vector in the documentation)
    for the data point being processed.
    Casting the coefficients into a matrix makes the calculation simple.
    */
    Matrix coef = state.coef;
    coef.resize(numCategories, state.widthOfX/numCategories);

    //Store the intermediate calculations because we'll reuse them in the LLH
    ColumnVector t1 = x; //t1 is vector of size state.widthOfX
    t1 = coef*x;
    /*
        Note: The above 2 lines could have been written as:
        ColumnVector t1 = -coef*x;

        but this creates warnings. These warnings are somehow related to the factor
        that x is an immutable type.
    */

    ColumnVector t2 = t1.array().exp();
    double t3 = 1 + t2.sum();
    ColumnVector pi = t2/t3;
    //The gradient matrix has numCatergories rows and widthOfX columns
    Matrix grad = -yVec * x.transpose() + pi * x.transpose();
    //We cast the gradient into a vector to make the math easier.
    grad.resize(state.widthOfX, 1);
    for (int i = 0; i < state.widthOfX; i++)
    {
		state.meat_half(0,i) += grad(i);
	}

	// Compute the 'a' matrix.
    Matrix a(numCategories,numCategories);
    Matrix piDiag = pi.asDiagonal();
    a = pi * pi.transpose() - piDiag;

    //Start the Hessian calculations
    //Matrix X_transp_AX(numCategories * state.widthOfX, numCategories * state.widthOfX);
 	Matrix X_transp_AX( (int)state.widthOfX, (int)state.widthOfX);
    /*