Drexel dragonThe Math ForumDonate to the Math Forum



Search All of the Math Forum:

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


Math Forum » Discussions » sci.math.* » sci.stat.math.independent

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

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   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
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

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



Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© Drexel University 1994-2014. All Rights Reserved.
The Math Forum is a research and educational enterprise of the Drexel University School of Education.