Function clustered_compute_stats

src/modules/regress/clustered_errors.cpp:184–262

Compute the stats from coef and errs

// Compute the stats from coef and errs
AnyType clustered_compute_stats (AnyType& args,
                                 void (*func)(
                                     MutableNativeColumnVector&,
                                     MutableNativeColumnVector&,
                                     int, int), bool ismlogr)
{
    //const MappedColumnVector& coef = args[0].getAs<MappedColumnVector>();
    ColumnVector coef;
    if (ismlogr){
	    MappedMatrix coefMat = args[0].getAs<MappedMatrix>();
        Matrix mat = coefMat;
        mat.transposeInPlace();
        mat.resize(coefMat.size(), 1);
        coef = mat;
    }else{
        coef = args[0].getAs<MappedColumnVector>();
    }
    const MappedColumnVector& meatvec = args[1].getAs<MappedColumnVector>();
    const MappedColumnVector& breadvec = args[2].getAs<MappedColumnVector>();
    int mcluster = args[3].getAs<int>();
    int numRows = args[4].getAs<int>();
    int k = static_cast<int>(coef.size());
    Matrix bread(k,k);
    Matrix meat(k,k);
    int count = 0;

    for (int i = 0; i < k; i++)
        for (int j = 0; j < k; j++)
        {
            meat(i,j) = meatvec(count);
            bread(i,j) = breadvec(count);
            count++;
        }

    if (mcluster == 1) {
        //throw std::domain_error ("Clustered variance error: Number of clusters cannot be smaller than 2!");
        warning("Clustered variance error: Number of clusters cannot be smaller than 2!");
        return Null();
    }
    double dfc = (mcluster / (mcluster - 1.)) * ((numRows - 1.) / (numRows - k));

    SymmetricPositiveDefiniteEigenDecomposition<Matrix> decomposition(
        bread, EigenvaluesOnly, ComputePseudoInverse);
    Matrix inverse_of_bread = decomposition.pseudoInverse();
    Matrix cov(k, k);
    cov = inverse_of_bread * meat * inverse_of_bread;

    MutableNativeColumnVector errs;
    MutableNativeColumnVector stats;
    MutableNativeColumnVector pValues;
    Allocator& allocator = defaultAllocator();

    errs.rebind(allocator.allocateArray<double>(k));
    stats.rebind(allocator.allocateArray<double>(k));
    pValues.rebind(allocator.allocateArray<double>(k));
    for (int i = 0; i < k; i++)
    {
        if (inverse_of_bread(i,i) < 0)