Newer
Older
/* glm.cc
Mark Woolrich, FMRIB Image Analysis Group
Copyright (C) 1999-2000 University of Oxford */
/* CCOPYRIGHT */
#include "glm.h"
#include "miscmaths.h"
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
void Glm::setData(Matrix& p_y, Matrix& p_x, Matrix& p_contrasts)
{
y = &p_y;
x = &p_x;
numTS = p_y.Ncols();
sizeTS = p_y.Nrows();
contrasts = &p_contrasts;
SetContrast(1);
}
const Matrix& Glm::ComputeResids()
{
Tracer ts("ComputeResids");
Matrix& d = *x;
Computeb();
// r = (I - x(x'x)-1x')y
Matrix I(sizeTS, sizeTS);
Identity(I);
r = (I-d*inv_xx*d.t())*(y->Column(1));
return r;
}
const float Glm::Computecb()
{
Tracer ts("Computecb");
cb = (c.t()*b).AsScalar();
return cb;
}
const ColumnVector& Glm::Computeb()
{
Tracer ts("Computeb");
Matrix& d = *x;
// inv_xx = inv(x'x)
inv_xx = (d.t()*d).i();
b = inv_xx*d.t()*(y->Column(1));
return b;
}
const float Glm::ComputeVar()
{
Tracer ts("ComputeVar");
// var = e*var_on_e
// e is the estimate of the variance of the timeseries, sigma^2
var = (r.Column(1).t()*r.Column(1)*c.t()*inv_xx*c).AsScalar()/sizeTS;
return var;
}
#ifndef NO_NAMESPACE
}
#endif