Skip to content

Commit

Permalink
Fix some wording in week3.qmd
Browse files Browse the repository at this point in the history
  • Loading branch information
teddygroves committed Mar 18, 2024
1 parent 3c08e1d commit d54d2e2
Show file tree
Hide file tree
Showing 4 changed files with 95 additions and 91 deletions.
44 changes: 22 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 @@ -367,7 +367,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([ 6., 45., 377., 1127., 1581., 1315., 503., 105., 9.,\n 2.]),\n array([0. , 0.43380605, 0.8676121 , 1.30141815, 1.7352242 ,\n 2.16903025, 2.6028363 , 3.03664235, 3.4704484 , 3.90425445,\n 4.3380605 ]),\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.",
"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., 7., 110., 545., 1309., 1525., 1033., 329., 55.,\n 4.]),\n array([-0.34331577, 0.09770966, 0.5387351 , 0.97976053, 1.42078597,\n 1.8618114 , 2.30283684, 2.74386228, 3.18488771, 3.62591315,\n 4.06693858]),\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 d54d2e2

Please sign in to comment.