Skip to content
Snippets Groups Projects
Commit f6087332 authored by Stephen Smith's avatar Stephen Smith
Browse files

change some .i() matrix inverses to use pinv() instead, in order to allow film...

change some .i() matrix inverses to use pinv() instead, in order to allow film to work with rank deficient design matrices
parent f5b330a9
No related branches found
No related tags found
No related merge requests found
...@@ -84,7 +84,8 @@ namespace FILM { ...@@ -84,7 +84,8 @@ namespace FILM {
int batch_pos = 1; int batch_pos = 1;
// pinv(x) = inv(x'x)x' // pinv(x) = inv(x'x)x'
pinv_x = (x.t()*x).i()*x.t(); //pinv_x = (x.t()*x).i()*x.t();
pinv_x = pinv(x);
// R = I - x*pinv(x) // R = I - x*pinv(x)
Matrix I(sizeTS, sizeTS); Matrix I(sizeTS, sizeTS);
...@@ -152,7 +153,7 @@ namespace FILM { ...@@ -152,7 +153,7 @@ namespace FILM {
Tracer ts("Glim::UseGlobalVrow"); Tracer ts("Glim::UseGlobalVrow");
// var/e = inv(x'x)x'*V*x*inv(x'x) // var/e = inv(x'x)x'*V*x*inv(x'x)
Matrix corr = pinv_x*V*x*(x.t()*x).i(); Matrix corr = pinv_x*V*pinv_x.t();
for(int i = 1; i <= numTS; i++) for(int i = 1; i <= numTS; i++)
{ {
......
...@@ -30,7 +30,8 @@ namespace FILM { ...@@ -30,7 +30,8 @@ namespace FILM {
void GlimGls::setData(const ColumnVector& y, const Matrix& x, const int ind) void GlimGls::setData(const ColumnVector& y, const Matrix& x, const int ind)
{ {
// compute b // compute b
Matrix inv_xx = (x.t()*x).i(); //Matrix inv_xx = (x.t()*x).i();
Matrix inv_xx = pinv(x.t()*x);
b.Column(ind) = inv_xx*x.t()*y; b.Column(ind) = inv_xx*x.t()*y;
// compute r // compute r
......
...@@ -54,7 +54,8 @@ namespace FILM { ...@@ -54,7 +54,8 @@ namespace FILM {
Matrix& d = *x; Matrix& d = *x;
// inv_xx = inv(x'x) // inv_xx = inv(x'x)
inv_xx = (d.t()*d).i(); // inv_xx = (d.t()*d).i();
inv_xx = pinv(d.t()*d);
b = inv_xx*d.t()*(y->Column(1)); b = inv_xx*d.t()*(y->Column(1));
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment