|
| 1 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
1 | 2 | % Particle Metropolis-Hastings (PMH) for the SV model |
| 3 | +% (c) Johan Dahlin 2017 under MIT license <liu@johandahlin.com.nospam> |
| 4 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
2 | 5 | function[theta, xHatFiltered] = particleMetropolisHastingsSVmodel(observations, initialParameters, noParticles, noIterations, stepSize) |
3 | 6 | noObservations = length(observations); |
4 | 7 |
|
|
14 | 17 | theta(1, :) = initialParameters; |
15 | 18 | [xHatFiltered(1, :), logLikelihood(1)] = particleFilterSVmodel(observations, theta(1, :), noParticles); |
16 | 19 |
|
17 | | - % Main loop |
18 | | - for k = 2:noIterations |
19 | | - |
| 20 | + for k = 2:noIterations |
20 | 21 | % Propose a new parameter |
21 | 22 | thetaProposed(k, :) = mvnrnd(theta(k-1, :), stepSize); |
22 | 23 |
|
|
25 | 26 | [xHatFilteredProposed(k, :), logLikelihoodProposed(k)] = particleFilterSVmodel(observations, thetaProposed(k, :), noParticles); |
26 | 27 | end |
27 | 28 |
|
28 | | - % Compute the acceptance probability |
| 29 | + % Compute the acceptance probability (reject if unstable) |
29 | 30 | prior = dnorm(thetaProposed(k, 1), 0, 1); |
30 | 31 | prior = prior - dnorm(theta(k - 1, 1), 0, 1); |
31 | 32 | prior = prior + dnorm(thetaProposed(k, 2), 0.95, 0.05); |
32 | 33 | prior = prior - dnorm(theta(k - 1, 2), 0.95, 0.05); |
33 | 34 | prior = prior + dgamma(thetaProposed(k, 3), 2, 10); |
34 | 35 | prior = prior - dgamma(theta(k - 1, 3), 2, 10); |
35 | | - |
36 | 36 | likelihoodDifference = logLikelihoodProposed(k) - logLikelihood(k - 1); |
37 | 37 | acceptProbability = exp(prior + likelihoodDifference); |
38 | | - |
39 | 38 | acceptProbability = acceptProbability * (abs(thetaProposed(k, 2)) < 1.0); |
40 | 39 | acceptProbability = acceptProbability * (thetaProposed(k, 3) > 0.0); |
41 | 40 |
|
42 | 41 | % Accept / reject step |
43 | 42 | uniformRandomVariable = unifrnd(0, 1); |
44 | | - |
45 | 43 | if (uniformRandomVariable < acceptProbability) |
46 | 44 | % Accept the parameter |
47 | 45 | theta(k, :) = thetaProposed(k, :); |
|
69 | 67 | end |
70 | 68 | end |
71 | 69 |
|
| 70 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
72 | 71 | % Helper for computing the logarithm of N(x; mu, sigma^2) |
| 72 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
73 | 73 | function[out] = dnorm(x, mu, sigma) |
74 | 74 | out = -0.5 .* log(2 * pi) - 0.5 .* log(sigma.^2) - 0.5 ./ sigma.^2 .* (x - mu).^2; |
75 | 75 | end |
76 | 76 |
|
| 77 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
77 | 78 | % Helper for computing the logarithm of Gamma(x; a, b) with mean a/b |
| 79 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
78 | 80 | function[out] = dgamma(x, a, b) |
79 | 81 | out = a * log(b) - gammaln(a) + (a-1) * log(x) - b * x; |
80 | 82 | end |
0 commit comments