0
$\begingroup$

I am trying to implement Metropolis–Hastings algorithm to find parameters. I basically followed this PDF. It seems it works but gives me off parameters. Note: This is the Metropolis algorithm not Metropolis–Hastings algorithm. I know there is MCMC package is available but I want to understand method. Can someone help me out to modified it? Thanks in advance.

Note: Cross-posted at http://community.wolfram.com/groups/-/m/t/1156904

a0 = 2; (*parameter 1*)
b0 = 5;  (*parameter 2*)
y[x_] := a0 x + b0 (*model*)

data = y@Range@10;

prior = PDF[NormalDistribution[0, 1], a] PDF[NormalDistribution[0, 1], b]
likelihood = Likelihood[NormalDistribution[a, b], data]


posterior[a_, b_] = prior *likelihood

ratio[{a1_, b1_}, {a2_, b2_}] =  posterior[a2, b2]/posterior[a1, b1] // Simplify


proposal[\[Theta]_,  std_] := \[Theta] + (RandomReal[NormalDistribution[0, #]] & /@ std)

test[\[Phi]_] := If[\[Phi][[2]] > 0, True, False]

metroPolisStep[\[Theta]_, std_] :=With[{\[Phi] = proposal[\[Theta], std], r = RandomReal[]},If[test[\[Phi]],  If[r <= ratio[\[Theta], \[Phi]], \[Phi], \[Theta]], \[Theta]]]

metroPolis[initialState_, std_, steps_] :=  NestList[metroPolisStep[#, std] &, initialState, steps]

numStep = 10000;

burnin = Ceiling[numStep* 10/100];

metro = metroPolis[{2, 1}, {0.5, 1}, numStep];

param1 = Drop[metro[[All, 1]], burnin];

param2 = Drop[metro[[All, 2]], burnin];

ListLinePlot[param1, PlotStyle -> PointSize[Tiny], AspectRatio -> 0.2,  ImageSize -> 400, PlotRange -> All]

ListLinePlot[param2, PlotStyle -> PointSize[Tiny], AspectRatio -> 0.2, ImageSize -> 400, PlotRange -> All]


{SmoothHistogram[param1, Axes -> {True, False}, ImageSize -> 300,    PlotLabel -> "a", PlotRange -> All],   SmoothHistogram[param2, Axes -> {True, False}, ImageSize -> 300,    PlotLabel -> "b", PlotRange -> All]} // GraphicsRow
$\endgroup$
  • $\begingroup$ Since SE does not support attachments, you could have at least edited your question here to link to the PDF instead of just doing copy-pasta. $\endgroup$ – J. M. will be back soon Aug 1 '17 at 5:32
  • $\begingroup$ Sorry, I did that. $\endgroup$ – OkkesDulgerci Aug 1 '17 at 15:11
  • 2
    $\begingroup$ I'm confused by your text. You want to implement the Metropolis-Hastings algorithm but you state "Note: This is the Metropolis algorithm not Metropolis–Hastings algorithm." $\endgroup$ – JimB Aug 1 '17 at 15:20
  • $\begingroup$ Main difference is Metropolis algorithm assumes dist. is symmetric and Metropolis-Hastings algorithm is more general case. $\endgroup$ – OkkesDulgerci Aug 1 '17 at 15:27
  • $\begingroup$ Would this help? 12000.org/my_notes/hastings_metropolis/index.htm Note there is a link for download of a notebook. $\endgroup$ – bill s Aug 1 '17 at 17:55
1
$\begingroup$

The following does not answer the OP's question directly, in that it does not provide modifications of the code presented. The purpose of this "answer" is to provide a clear statement of the Metropolis-Hastings algorithm and its relation to the Metropolis algorithm in hopes that this would aid the OP in modifying the code him- or herself. My impression from reading the OP's question is that there may be a misunderstanding on the OP's part, but perhaps I am wrong and this "answer" should be removed.

Given data (vector) $y$ and unknown parameter (vector) $\theta$, the posterior distribution is given by $$ p(\theta|y) = \frac{p(y|\theta)\,p(\theta)}{p(y)} , $$ where $p(y|\theta)$ is the likelihood, $p(\theta)$ is the prior, and $p(y)$ is (often called) the marginal likelihood.

The Metropolis-Hastings algorithm can be summarized as follows. Let $q(\theta'|\theta)$ denote the density for the proposal distribution for $\theta'$ conditioned on $\theta$. Let $\theta^{(r)}$ denote the current state of the MCMC chain. Then \begin{equation} \theta^{(r+1)} = \begin{cases} \theta' & \mathcal{R} \ge u \\ \theta^{(r)} & \text{otherwise} \end{cases} , \end{equation} where $u \sim \textsf{Uniform}(0,1)$ and \begin{equation} \mathcal{R} = \frac{p(\theta'|y)}{p(\theta^{(r)}|y)} \times \frac{q(\theta^{(r)}|\theta')}{q(\theta'|\theta^{(r)})} . \end{equation} Note that $p(y)$ cancels out of the first factor on the right-hand side. Also note that if \begin{equation} q(\theta|\theta') = q(\theta'|\theta) \qquad\text{for all $(\theta,\theta')$} \end{equation} then the second factor disappears and we have the Metropolis algorithm.

$\endgroup$
  • 1
    $\begingroup$ This doesn't seem to have anything to do with Mathematica. The OP wanted help with modifying the algorithm. A Mathematica implementation of what you wrote would be more in line with the question. $\endgroup$ – b3m2a1 Aug 14 '17 at 18:27
  • $\begingroup$ @mef Thanks for your post. I have read about theory but I need help for implementation. $\endgroup$ – OkkesDulgerci Aug 16 '17 at 3:43
0
$\begingroup$

I tried this but I am still not sure if I am doing right.

ClearAll["Global`*"]

SeedRandom[123]
m0 = 3;
b0 = 5;
std = 5;
sampleSize = 31;
y[x_] := m0 x + b0

data = Thread[{Range[0, 30], 
   y@Range[0, 30] + 
    RandomVariate[NormalDistribution[0, std], sampleSize]}]

Show[Plot[y[x], {x, 0, 30}, Frame -> True, Axes -> False], 
 ListPlot[data, PlotStyle -> Red]]

model[x_] := m x + b

prior = PDF[NormalDistribution[3, 1], m] PDF[NormalDistribution[5, 1],
    b]

Variance[data[[All, 2]]]

n = 31;
\[Sigma] = StandardDeviation[data[[All, 2]]]

likelihood = 
 Product[(1/Sqrt[2 \[Pi] (\[Sigma]^2) ]
     Exp[-(1/(
       2 \[Sigma]^2)) (data[[i, 2]] - model[data[[i, 1]]])^2]), {i, n}]

posterior[m_, b_] = prior likelihood

pdf = Compile[{{m, _Real}, {b, _Real}}, posterior[m, b], 
   Parallelization -> True];


mcmcfunComp = 
  Compile[{{s, _Real, 1}}, 
   Module[{xm, ym, proposal, xp, yp, p2, p1, proposalSigma = 0.2},
    {xm, ym} = s;
    xp = RandomReal[NormalDistribution[xm, proposalSigma]];
    yp = RandomReal[NormalDistribution[ym, proposalSigma]];
    p2 = pdf[xp, yp];
    p1 = pdf[xm, ym];
    proposal = p2/p1;
    If[RandomReal[] <= proposal, {xp, yp}, {xm, ym}]], 
   CompilationOptions -> {"InlineExternalDefinitions" -> True}, 
   Parallelization -> True];

sim = Drop[NestList[mcmcfunComp[#] &, {2, 4}, 100000], 10000];

ListLinePlot[#, PlotRange -> All, PlotStyle -> Thin, 
    AspectRatio -> 0.3, Frame -> True, 
    ImageSize -> Medium] & /@ {sim[[All, 1]], 
   sim[[All, 2]]} // GraphicsRow

{Histogram[sim[[All, 1]], Axes -> {True, False}, ImageSize -> 300, 
   PlotLabel -> "a", PlotRange -> All, 
   ChartElementFunction -> "FadingRectangle",
   ChartStyle -> Orange], 
  Histogram[sim[[All, 1]], Axes -> {True, False}, ImageSize -> 300, 
   PlotLabel -> "b", PlotRange -> All, 
   ChartElementFunction -> "FadingRectangle",
   ChartStyle -> Orange]} // GraphicsRow

{Mean@sim[[All, 1]], Mean@sim[[All, 2]]}

{SmoothHistogram[sim[[All, 1]], Axes -> {True, False}, 
   ImageSize -> 300, PlotLabel -> "m", PlotRange -> All], 
  SmoothHistogram[sim[[All, 2]], Axes -> {True, False}, 
   ImageSize -> 300, PlotLabel -> "b", 
   PlotRange -> All]} // GraphicsRow

Labeled[SmoothDensityHistogram[Thread[{sim[[All, 1]], sim[[All, 2]]}],
   "Oversmooth", "PDF", ColorFunction -> "ThermometerColors", 
  ImageSize -> 300, Mesh -> 3, MeshStyle -> Directive[White, Dashed], 
  GridLines -> Automatic, 
  PlotRange -> All], {Style["m ", 20, FontFamily -> "Times New Roman",
    Bold], Style["b ", 20, FontFamily -> "Times New Roman", 
   Bold]}, {Bottom, Left}]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.