pf_bootstrap.jl 4.85 KB
Newer Older
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
1
2
# module Tmp
using Devectorize
3
using StatsBase
4
function pf_bootstrap!(xp,w, y, N, g_density, f)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
5
    T = length(y)
6
7
8
    xp[:,1] = 2*σw*randn(N)
    wT = slice(w,:,1)
    xT = slice(xp,:,1)
9
    fill!(wT,-log(N))
10
11
    g(y[1],xT, wT)
    wT -= logsumexp(wT)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
12
13
14
    j = Array(Int64,N)
    #     tw = Float64(N)
    for t = 2:T
15
16
17
18
        xT1 = xT
        xT = slice(xp,:,t)
        wT1 = wT
        wT = slice(w,:,t)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
19
20
        # Resample
        if t%2 == 0
21
22
            resample_systematic!(j,wT1,N)
            f(xT,xT1[j],t-1)
23
            fill!(wT,-log(N))
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
24
        else # Resample not needed
25
26
            f(xT,xT1,t-1)
            copy!(wT,wT1)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
27
28
29
        end

        # Time update
30
        g(y[t],xT, wT)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
31
32

        # Normalize weights
33
        wT -= logsumexp(wT)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
34
35
36
37
38
39
        #         tw =  1/sum(w(:,t).*2)
    end

    return nothing
end

40
function pf_bootstrap(y, N, g_density, f)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
41
42
43
    T = length(y)
    xp = Array(Float64,(N,T))
    w = Array(Float64,(N,T))
44
45
46
    xp[:,1] = 2*σw*randn(N)
    wT = slice(w,:,1)
    xT = slice(xp,:,1)
47
    fill!(wT,-log(N))
48
49
    g(y[1],xT, wT)
    wT -= logsumexp(wT)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
50
51
52
    j = Array(Int64,N)
    #     tw = Float64(N)
    for t = 2:T
53
54
55
56
        xT1 = xT
        xT = slice(xp,:,t)
        wT1 = wT
        wT = slice(w,:,t)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
57
58
        # Resample
        if t%2 == 0
59
60
            resample_systematic!(j,wT1,N)
            f(xT,xT1[j],t-1)
61
            fill!(wT,-log(N))
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
62
        else # Resample not needed
63
64
            f(xT,xT1,t-1)
            copy!(wT,wT1)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
65
66
67
        end

        # Time update
68
        g(y[t],xT, wT)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
69
70

        # Normalize weights
71
        wT -= logsumexp(wT)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
72
73
74
75
76
77
78
79
80
        #         tw =  1/sum(w(:,t).*2)
    end

    return xp,w
end




81
function pf_bootstrap_nn(xp, w, wnn, y, N, g_density, f)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
82
    T = length(y)
83
84
85
    xp[:,1] = 2*σw*randn(N)
    wT = slice(w,:,1)
    xT = slice(xp,:,1)
86
    fill!(wT,-log(N))
87
88
89
    g(y[1],xT, wT)
    wnn[:,t] = wT
    wT -= logsumexp(wT)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
90
91
92
    j = Array(Int64,N)
    #     tw = Float64(N)
    for t = 2:T
93
94
95
96
        xT1 = xT
        xT = slice(xp,:,t)
        wT1 = wT
        wT = slice(w,:,t)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
97
98
        # Resample
        if t%2 == 0
99
100
            resample_systematic!(j,wT1,N)
            f(xT,xT1[j],t-1)
101
            fill!(wT,-log(N))
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
102
        else # Resample not needed
103
104
            f(xT,xT1,t-1)
            copy!(wT,wT1)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
105
106
107
108
109
110
111
        end

        # Time update
        g_density(y[t],xp[:,t], wT)
        wnn[:,t] = wT
        # Normalize weights
        #some strange error with this normalization method, I totally have no clue why it's working elsewhere but not here
112
113
        wT -= logsumexp(wT)
        #         @assert isapprox(nc,log(sum(exp(wnn[:,t]))))
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
114
115
116
117
118
119
120
        #         tw =  1/sum(w(:,t).*2)
    end

    return xp,w, wnn
end


121
function pf_aux_nn(xp, w, wnn, y, N, g_density, f)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
122
    T = length(y)
123
    xp[:,1] = 2*σw*randn(N)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
124
    w0 = fill(log(1/N),N)
125
    xp[:,1] = f(zeros(N),1)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
    wnn[:,1] = w0 + g_density(y[1],xp[:,1])
    @devec w[:,1] = wnn[:,1] - log(sum(exp(wnn[:,1])))
    j = Array(Int64,N)
    #     tw = Float64(N)
    for t = 2:T

        xpT = xp[:,t-1]
        f(xpT,t-1)

        wT = copy(w[:,t-1])
        lambda = g_density(y[t],xpT)
        wT += lambda
        wTe = exp(wT)
        wTe /= sum(wTe)

        # Resample
        if t%2 == 0
            resample_systematic_exp!(j,wTe)
            xpT = xp[j,t-1]
            wT = copy(w0)
            lambdaT = lambda[j]
        else # Resample not needed
            xpT = xp[:,t-1]
            #             wT = copy(w[:,t-1])
            lambdaT = lambda
        end

        # Time update
        f(xpT,t-1)
155
        xp[:,t] = xpT
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
        wT += g_density(y[t],xp[:,t])
        wnn[:,t] = wT - lambdaT
        # Normalize weights
        #some strange error with this normalization method, I totally have no clue why it's working elsewhere but not here
        # #         offset = maximum(wT)
        # #         normConstant = 0.0
        # #         for i = 1:N
        # #             normConstant += exp(wT[i]-offset)
        # #         end
        # #         w[:,t] = wT - log(normConstant)+offset
        w[:,t] = wnn[:,t] - log(sum(exp(wnn[:,t])))
        #         tw =  1/sum(w(:,t).*2)
    end

    return xp,w, wnn
end
172
173
174


function pf_CSMC!(xp,w, y, N, g_density, f, xc )
175
    xp[:,1] = 2*σw*randn(N)
176
    xp[N,1] = xc[1]
177
178
179
    T = length(y)
    wT = slice(w,:,1)
    xT = slice(xp,:,1)
180
    fill!(wT,-log(N))
181
182
    g(y[1],xT, wT)
    wT -= logsumexp(wT)
183
184
185
    j = Array(Int64,N)
    #     tw = Float64(N)
    for t = 2:T
186
187
188
189
        xT1 = xT
        xT = slice(xp,:,t)
        wT1 = wT
        wT = slice(w,:,t)
190
191
        # Resample
        if t%2 == 0
192
193
            resample_systematic!(j,wT1,N)
            f(xT,xT1[j],t-1)
194
            fill!(wT,-log(N))
195
        else # Resample not needed
196
197
            f(xT,xT1,t-1)
            copy!(wT,wT1)
198
199
        end
        xp[N,t] = xc[t]
200
201
        # Time update
        g(y[t],xT, wT)
202
203

        # Normalize weights
204
        wT -= logsumexp(wT)
205
206
207
208
209
        #         tw =  1/sum(w(:,t).*2)
    end

    return nothing
end