Math Forum Discussions

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

Plain Text

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);