-
Notifications
You must be signed in to change notification settings - Fork 1
/
pca.jl
70 lines (56 loc) · 1.35 KB
/
pca.jl
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
####
#
# EXAMPLE
#
####
using DyadicMatrices
using Winston
DM = DyadicMatrices
DV = DyadicVector
const L = 8
const K = 0
#srand(4)
alpha = 0.5
a = 5/2
b = 5/2 #hyper prior on si0
T = 1.
N = 500
# f(t) = 2sin(2pi ./(0.9t +.1)) #- 4hat(2t-1)
f(t) = 2sqrt(sin(4pi*t)+1) - 2t
tt = linspace(0, T, N)
sitrue = 2.
X = f(tt) + sitrue*randn(N)
n = len(L, K)
si0 = 4
si = std(diff(X)/sqrt(2)) # simple estimate for empirical bayes
mu = DM.mureg(tt, X, si, L, K)
mubar = DM.fetransf(f, 0.0, T, 8, K)
Sigma = DM.Sigmareg(tt, si, L, K)
mu = DM.pickupmu!(mu)
Sigma = pickupSigma(Sigma, K)
W = addprecision(Sigma, si0, alpha)
R = root(W)
theta0 = R'\(R\mu)
x = DM.drop(theta0)
println("Res: ", round(std(X.-evalfe(tt,x)),2))
println("si0 ", round(si0,2), " si ", round(si,2))
S = inv!(copy(R))
band = 1.64*sqrt(droppedvar(S))
thetar = R'\(R\mu + DV(randn(n),K))
xr = DM.drop(thetar)
n2 = n + 2 - 2K
if K == 0
dp = T/n:T/(n+1):T*(n-1)/n
else
dp = linspace(0, T, n)
end
if true
display(plot(tt, X, ".",
linspace(0, 1, length(mubar) ), mubar, "r",
dp, x.x, "g-",
dp, xr.x, "g.",
dp, x.+band, "b--", dp, x.-band, "b--",
yrange=[-6,6], linewidth=1))
print("Ok")
end
# readline()