This is a sample of a gaussian density using the "native" method:
μ
1x2 1 1
cov
2x2 1.000 0.800 0.800 1.000
disp 3 (asRow $ fst okt)
1x2 1.038 1.034
disp 3 (snd okt)
2x2 0.815 0.660 0.660 0.876
metropolis' :: Seed -> Int -> Int -> ℝ -> V -> (V -> ℝ) -> [MCMC]
metropolis' seed step burn σ x0 lprob = g $ scanl f mc0 deltas
where
f (MCMC x lp _) (δ,r)
| accept = MCMC x' lp' True
| otherwise = MCMC x lp False
where
x' = x+δ
lp' = lprob x'
ratio = lp' - lp
accept = ratio > 0 || log r < ratio
deltas = zip (infn seed (size x0) σ) (infu seed)
g = drop burn . map head . chunksOf step
mc0 = MCMC x0 (lprob x0) True
(deltas is an infinite list of tuples of gaussian perturbations of the appropriate
dimensions and uniform values)
The initial raw samples (including the "burn-in" time and without any thining):
The resulting "good" sample:
disp 3 (asRow $ fst ok)
1x2 1.341 1.180
disp 3 (snd ok)
2x2 1.272 1.017 1.017 1.288
(Convergence to the true parameters can be confirmed with a larger sample size.)
(test2,ar2) = info $ take 300 $ metropolis' seed step burn σ x0 logProb
where
seed = 888
step = 10
burn = 1000
σ = 1
x0 = vector [-3,3]
logProb = lgauss 2 (0.3^2) . scalar . norm_2
The acceptance ratio is 0.35 .
This is the example in the excellent post frequentism vs bayesianism.
The 5-3 advantage gives a maximum likelihood estimate $\hat p = 3/8= 0.375 $, and then a probability of three consecutive wins $\hat p^3$ = 0.05.
The full bayesian analisis gives a different result. For a uniform prior on $p$ the posterior is:
$$P(p|3-5)=\frac{p^3(1-p)^5}{B(4,6)}$$
The probability of tree consecutive wins with this distribution is
$$P(6-5) = \frac{B(7,6)}{B(4,6)} \simeq 0.091$$
We can try to replicate this result using Metropolis' MCMC. First we compute samples from the posterior distribution of $p$. This is a very simple one-dimensional case.
mapM_ (printf "%5.2f\n") $ take 10 $ map ((!0)) test4
0.36 0.28 0.36 0.12 0.24 0.28 0.77 0.35 0.46 0.47
The acceptance ratio is 0.202 . The histogram is consistent with the analytic solution derived above:
And then we compute the empirical expectation of the probability of three more wins:
printf "%.3f" $ average (map ((**3).(!0)) test4)
0.094