Skip to content

Commit

Permalink
updating accept_reject value
Browse files Browse the repository at this point in the history
  • Loading branch information
NicholasCowie committed Mar 5, 2024
1 parent a4931cb commit 90c5b38
Show file tree
Hide file tree
Showing 4 changed files with 85 additions and 76 deletions.
46 changes: 24 additions & 22 deletions docs/metropolis-hastings.html

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion docs/search.json
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,7 @@
"href": "metropolis-hastings.html#coding-the-metropolis-hastings-in-practice",
"title": "Metropolis Hastings",
"section": "Coding the Metropolis Hastings in practice",
"text": "Coding the Metropolis Hastings in practice\n\nPart 1: sampling from a normal distribution\nGiven a \\(p(x) = N(2, 0.5)\\) how would draw samples using the Metropolis-Hastings algorithm?\n\nChoose proposal value\nEvaluate probability ratio\nAccept-Reject\nincrement step\n\n\nDefine probability density function\n\nimport numpy as np\nfrom scipy.stats import norm\n\ndef prob(x):\n return norm.pdf(x, 2, 0.5)\n\n\n\nDefine proposal distribution\n\ndef proposal(x):\n return norm.rvs(x, 1, 1)[0]\n\n\n\nInitialise sampler\n\ncurrent = 0.0\nsamples = [current]\n\n\n\nSample from distribution\n\nfor i in range(10000):\n prop = proposal(current)\n accept_reject = prob(prop)/prob(current)\n if accept_reject > 1:\n samples.append(prop)\n current = prop\n else:\n cutoff = np.random.rand(1)[0]\n if accept_reject > cutoff:\n samples.append(prop)\n current = prop\n\n\n\nPlot distribution\n\nimport matplotlib.pyplot as plt\nplt.hist(samples)\n\n(array([ 4., 23., 158., 627., 1159., 1379., 1029., 476., 113.,\n 17.]),\n array([-0.03049052, 0.34894132, 0.72837315, 1.10780499, 1.48723683,\n 1.86666866, 2.2461005 , 2.62553234, 3.00496418, 3.38439601,\n 3.76382785]),\n <BarContainer object of 10 artists>)\n\n\n\n\n\n\n\n\n\n\n\nTrace plot\n\ndraw = [draw for draw, _ in enumerate(samples)]\nplt.plot(draw, samples)\n\n\n\n\n\n\n\n\n\n\n\nPart 2: determining mean and standard deviation from data\nI suggest using logs due to numerical issues. Here’s an example function which you can use to evaluate the probability of the data.\n\ndef eval_prob(data, mu, sigma):\n return np.log(norm.pdf(data,mu,sigma)).sum()\n\nHere’s also a multivariate random number generator to generate proposals.\n\ndef proposal_multi(mu, sigma):\n mean = [mu, sigma]\n cov = [[0.2, 0], [0, 0.2]] # diagonal covariance\n return np.random.multivariate_normal(mean, cov, 1)[0]\n\nHere is how you’d call the proposal function\nprop_mu, prop_sigma = proposal_multi(current_mu, current_sigma)\nYou should sample a 95% interval including a \\(\\mu = 5\\) and a \\(\\sigma = 0.2\\). This may be difficult at first to sample and I would recommend initialising at these values.",
"text": "Coding the Metropolis Hastings in practice\n\nPart 1: sampling from a normal distribution\nGiven a \\(p(x) = N(2, 0.5)\\) how would draw samples using the Metropolis-Hastings algorithm?\n\nChoose proposal value\nEvaluate probability ratio\nAccept-Reject\nincrement step\n\n\nDefine probability density function\n\nimport numpy as np\nfrom scipy.stats import norm\n\ndef prob(x):\n return norm.pdf(x, 2, 0.5)\n\n\n\nDefine proposal distribution\n\ndef proposal(x):\n return norm.rvs(x, 1, 1)[0]\n\n\n\nInitialise sampler\n\ncurrent = 0.0\nsamples = [current]\n\n\n\nSample from distribution\n\nfor i in range(10000):\n prop = proposal(current)\n accept_reject = prob(prop)/prob(current)\n if accept_reject > 1:\n samples.append(prop)\n current = prop\n else:\n cutoff = np.random.rand(1)[0]\n if accept_reject > cutoff:\n samples.append(prop)\n current = prop\n\n\n\nPlot distribution\n\nimport matplotlib.pyplot as plt\nplt.hist(samples)\n\n(array([ 6., 30., 204., 606., 1071., 1257., 1130., 536., 161.,\n 13.]),\n array([0. , 0.36895147, 0.73790294, 1.10685441, 1.47580588,\n 1.84475735, 2.21370882, 2.58266029, 2.95161176, 3.32056323,\n 3.6895147 ]),\n <BarContainer object of 10 artists>)\n\n\n\n\n\n\n\n\n\n\n\nTrace plot\n\ndraw = [draw for draw, _ in enumerate(samples)]\nplt.plot(draw, samples)\n\n\n\n\n\n\n\n\n\n\n\nPart 2: determining mean and standard deviation from data\nI suggest using logs due to numerical issues. Here’s an example function which you can use to evaluate the probability of the data.\n\ndef eval_prob(data, mu, sigma):\n return np.log(norm.pdf(data,mu,sigma)).sum()\n\nHere’s also a multivariate random number generator to generate proposals.\n\ndef proposal_multi(mu, sigma):\n mean = [mu, sigma]\n cov = [[0.2, 0], [0, 0.2]] # diagonal covariance\n return np.random.multivariate_normal(mean, cov, 1)[0]\n\nHere is how you’d call the proposal function\nprop_mu, prop_sigma = proposal_multi(current_mu, current_sigma)\nthe accept_reject probability should also be updated to account for the log-probability\naccept_reject = np.exp(prob_prop - prob_current)\nYou should sample a 95% interval including a \\(\\mu = 5\\) and a \\(\\sigma = 0.2\\). This may be difficult at first to sample and I would recommend initialising at these values.",
"crumbs": [
"Course materials",
"Metropolis Hastings"
Expand Down
Loading

0 comments on commit 90c5b38

Please sign in to comment.