-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathpc_regession.go
93 lines (75 loc) · 1.71 KB
/
pc_regession.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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
87
88
89
90
91
92
93
package glasso
/*
import (
"github.com/gonum/matrix/mat64"
)
type PCR struct {
x *DataFrame
z *DataFrame
m int
}
func NewPCR(x *DataFrame, m int) *PCR {
z := mat64.DenseCopyOf(x.data)
r, c := z.Dims()
// need to standardize x for best results
x.Standardize()
return &PCR{
x: x,
z: &DataFrame{z, r, c, nil},
m: m,
}
}
/*
func (p *PCR) GetPrinComps(k float64) *mat64.Dense {
epsilon := math.Pow(2, -52.0)
small := math.Pow(2, -966.0)
svd := mat64.SVD(mat64.DenseCopyOf(p.x.data), epsilon, small, false, true)
// Sigma is a square + diagonal matrix
eigenvalues := svd.Sigma
l := len(S)
sum := 0.0
for i := 0; i < r; i++ {
sum += eigenvalues[i]
}
for i := 0; i < r; i++ {
eigenvalues[i] /= sum
}
varexplained := 0.0
r, c := p.x.rows, p.x.cols
if k > c {
log.Println("cannot have k (cutoff) > # columns")
return nil
}
// projection matrix W:
W := mat64.NewDense(c, k, nil)
// eigenvectors
V := svd.V
for i := 0; i < k; i++ {
var col []float64
V.Col(col, i)
W.SetCol(i, col)
varexplained += eigenvalues[i]
}
log.Printf(`Projection matrix created. \%%d of data explained`, varexplained)
out := &mat64.Dense{}
out.Mul(p.z.data, W)
return out
}
// Principal component regression forms the derived input columns zm = Xvm
// and then regresses y on z1, z2, . . . , zM for some M <= p. Since the zm
// are orthogonal, this regression is just a sum of univariate regressions:
//
// y_hat = y_bar + sum \theta_m z_m
// where \theta_m = <z_m, y> / <z_m, z_m>
//
func (p *PCR) Train(y []float64) error {
//ybar := mean(y)
theta := make([]float64, p.m)
for i := 0; i < p.m; i++ {
z_m := p.z.data.Col(nil, i)
theta[i] = sum(prod(z_m, y)) / sum(prod(z_m, z_m))
}
// ...
return nil
}
*/