Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Topic: Perl code for taking inv[X'X] from Ivo Welch's Statistics::Regression
to b,u,v,C via Jonathan Leto's Math::MatrixReals

Replies: 1   Last Post: Aug 28, 2012 1:25 PM

 Messages: [ Previous | Next ]
 Halitsky Posts: 600 Registered: 2/3/09
Perl code for taking inv[X'X] from Ivo Welch's Statistics::Regression
to b,u,v,C via Jonathan Leto's Math::MatrixReals

Posted: Aug 28, 2012 1:01 PM

Below is some dirty, on the fly, BUT tested and working code for
taking the inv[X'X] matrix returned by Ivo Welch's Perl
Statistics::Regression module and obtaining b, u, v, and C using
Jonathan Leto's PERL module Math:MatrixReals. (Following Ray
Koopman's notation, "b" = AX'y = OLS estimate of B, u = y-xb = sample
residuals, v = u'u/(n-k) = sample estimate of V, C = vA (estimate of
covariance matrix of B.)

Two things are important to note.

First, Ivo returns the ARRAY REFERENCE to inv[X'X] as the (zero-based)
(k+1)th entry in the array returned by his standarderrors method for a
multiple linear regression involving k predictors. For example, if
you ask Ivo's module to regress 1 DV on two predictors, the array
returned by this method will be:

\$array[0] = intercept
\$array[1] = slope of first predictor
\$array[2] = slope of second predictor
\$array[3] = array reference to inv[X'X] matrix

Second, the form of the inv[X'X] matrix in Ivo's module is NOT
compatible with the form of matrices demanded by Jonathan's module.
Hence, the code below does a little dance to convert one format to the
other.

Code:

use Statistics::Regression;
use Math::MatrixReal;
my \$reg = Statistics::Regression->new( "Leu", [ "const", "lncove",
"lncovu" ] );
my %d;
for (\$Ncntr=0;\$Ncntr<\$obsN;\$Ncntr++)
{
# send input rows to Ivo's module
%d = ();
\$dL = \$Ls[\$Ncntr];
\$d{const} = 1.0;
\$d{lncove} = \$es[\$Ncntr];
\$d{lncovu} = \$us[\$Ncntr];
\$d{ignored} = "anything else";
\$reg->include( \$dL, \%d );

#construct string for matrix X for use later after call to Ivo's
module
\$stringX = \$stringX . '[ ' . \$d{const} . ' ' .
\$d{lncove} . ' ' .
\$d{lncovu} . ' ]' . "\n";

#construct string for matrix Y for use later after call to Ivo's
module
\$stringY = \$stringY . '[ ' . \$dL . ' ]' . "\n";

}

\$mX = Math::MatrixReal->new_from_string(\$stringX);
\$mY = Math::MatrixReal->new_from_string(\$stringY);
\$mXp = new Math::MatrixReal(3,\$obsN);
\$mXp->transpose(\$mX);

my @theta = \$reg->theta();
my @se;
my @xpxinv;
my \$sigmasq;
@se = \$reg->standarderrors();

\$xpxinvref = \$se[3];

for (my \$i=1; \$i<=3; ++\$i)
{
\$stringA = \$stringA . '[ ';
for (my \$j=1; \$j<=3; ++\$j)
{
\$stringA = \$stringA . \$\$xpxinvref[ui(\$i,\$j)] . ' ';
}
\$stringA = \$stringA . ']' . "\n";
}

\$mA = Math::MatrixReal->new_from_string(\$stringA);

\$mXpy = \$mXp->multiply(\$mY);

\$mb = \$mA->multiply(\$mXpy);

\$mXb = \$mX->multiply(\$mb);

\$mu = new Math::MatrixReal(\$obsN,1);
\$mu->subtract(\$mY,\$mXb);

\$mup = new Math::MatrixReal(1,\$obsN);
\$mup->transpose(\$mu);

\$mupu = \$mup->multiply(\$mu);

\$df = 1 / ( \$obsN - 3 );
\$mv = new Math::MatrixReal(1,1);
\$mv->multiply_scalar(\$mupu,\$df);

\$v = \$mv->element(1,1);
\$mC = new Math::MatrixReal(3,3);
\$mC->multiply_scalar(\$mA,\$v);

Date Subject Author
8/28/12 Halitsky
8/28/12 Halitsky