From 90c5b389b5f8f40ed455d64546a2e1df421bd11c Mon Sep 17 00:00:00 2001 From: NicholasCowie Date: Tue, 5 Mar 2024 12:55:59 +0100 Subject: [PATCH] updating accept_reject value --- docs/metropolis-hastings.html | 46 ++++++------- docs/search.json | 2 +- docs/week3.html | 106 +++++++++++++++--------------- materials/metropolis-hastings.qmd | 7 ++ 4 files changed, 85 insertions(+), 76 deletions(-) diff --git a/docs/metropolis-hastings.html b/docs/metropolis-hastings.html index a094742..9b9595e 100644 --- a/docs/metropolis-hastings.html +++ b/docs/metropolis-hastings.html @@ -5007,7 +5007,7 @@

Markov Chains

For example:

The transitions can be measured as discrete time steps with the following matrix representation

-
+
import numpy as np
 T = np.matrix([[0.1, 0.1, 0.8], [0.5, 0.1, 0.4], [0.5, 0.2, 0.3]])
 print(T)
@@ -5018,7 +5018,7 @@

Markov Chains

Given an initial starting position

-
+
v0 = np.matrix([0.1, 0.4, 0.5])
 print(v0)
@@ -5026,7 +5026,7 @@

Markov Chains

We can simulate the probabilities of the next step given the transition matrix.

-
+
v1 = v0*T
 print(v1)
@@ -5034,7 +5034,7 @@

Markov Chains

Following this again we can simulate the states after two steps

-
+
v2 = v1*T
 print(v2)
@@ -5042,21 +5042,21 @@

Markov Chains

There’s a convenient way to calculate the next step given the starting condition.

-
+
print(v0*T**2)
[[0.316 0.139 0.545]]

What happens if we continue doing this for a long time?

-
+
print(v0*T**100)
[[0.35714286 0.14935065 0.49350649]]

And how does this change when taking the next step?

-
+
print(v0*T**101)
[[0.35714286 0.14935065 0.49350649]]
@@ -5121,7 +5121,7 @@

Define probability density function

-
+
import numpy as np
 from scipy.stats import norm
 
@@ -5131,21 +5131,21 @@ 
Define

Define proposal distribution

-
+
def proposal(x):
   return norm.rvs(x, 1, 1)[0]

Initialise sampler

-
+
current = 0.0
 samples = [current]

Sample from distribution

-
+
for i in range(10000):
     prop = proposal(current)
     accept_reject = prob(prop)/prob(current)
@@ -5161,21 +5161,21 @@ 
Sample from distr

Plot distribution

-
+
import matplotlib.pyplot as plt
 plt.hist(samples)
-
(array([   4.,   23.,  158.,  627., 1159., 1379., 1029.,  476.,  113.,
-          17.]),
- array([-0.03049052,  0.34894132,  0.72837315,  1.10780499,  1.48723683,
-         1.86666866,  2.2461005 ,  2.62553234,  3.00496418,  3.38439601,
-         3.76382785]),
+
(array([   6.,   30.,  204.,  606., 1071., 1257., 1130.,  536.,  161.,
+          13.]),
+ array([0.        , 0.36895147, 0.73790294, 1.10685441, 1.47580588,
+        1.84475735, 2.21370882, 2.58266029, 2.95161176, 3.32056323,
+        3.6895147 ]),
  <BarContainer object of 10 artists>)
-

+

@@ -5183,13 +5183,13 @@

Plot distribution

Trace plot

-
+
draw = [draw for draw, _ in enumerate(samples)]
 plt.plot(draw, samples)
-

+

@@ -5199,12 +5199,12 @@

Trace plot

Part 2: determining mean and standard deviation from data

I suggest using logs due to numerical issues. Here’s an example function which you can use to evaluate the probability of the data.

-
+
def eval_prob(data, mu, sigma):
     return np.log(norm.pdf(data,mu,sigma)).sum()

Here’s also a multivariate random number generator to generate proposals.

-
+
def proposal_multi(mu, sigma):
     mean = [mu, sigma]
     cov = [[0.2, 0], [0, 0.2]]  # diagonal covariance
@@ -5212,6 +5212,8 @@ 
\(\mu = 5\) and a \(\sigma = 0.2\). This may be difficult at first to sample and I would recommend initialising at these values.
diff --git a/docs/search.json b/docs/search.json index 6effb14..024d198 100644 --- a/docs/search.json +++ b/docs/search.json @@ -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" diff --git a/docs/week3.html b/docs/week3.html index fc0ba73..9333c64 100644 --- a/docs/week3.html +++ b/docs/week3.html @@ -5132,7 +5132,7 @@

Results

Example

We’ll go through some diagnostics using arviz.

-
+
import arviz as az
 import numpy as np
 import xarray as xr
@@ -5148,8 +5148,8 @@ 
Example

               
    
               
             
  • 
-                  
-                  
+                  
+                  
                   
    
                   
    
                       
      
@@ -5476,7 +5476,7 @@ 
      Example
      
     inference_library:          pymc3
     inference_library_version:  3.9.2
     sampling_time:              18.096983432769775
-    tuning_steps:               1000
  • created_at :
    2020-07-24T18:15:12.191355
    arviz_version :
    0.9.0
    inference_library :
    pymc3
    inference_library_version :
    3.9.2
    sampling_time :
    18.096983432769775
    tuning_steps :
    1000

  • - - + +
      @@ -5833,20 +5833,20 @@

      Example

      created_at: 2020-07-24T18:15:12.449843 arviz_version: 0.9.0 inference_library: pymc3 - inference_library_version: 3.9.2
  • created_at :
    2020-07-24T18:15:12.449843
    arviz_version :
    0.9.0
    inference_library :
    pymc3
    inference_library_version :
    3.9.2

  • - - + +
      @@ -6164,20 +6164,20 @@

      Example

      created_at: 2020-07-24T18:15:12.448264 arviz_version: 0.9.0 inference_library: pymc3 - inference_library_version: 3.9.2
  • created_at :
    2020-07-24T18:15:12.448264
    arviz_version :
    0.9.0
    inference_library :
    pymc3
    inference_library_version :
    3.9.2

  • - - + +
      @@ -6505,17 +6505,17 @@

      Example

      inference_library: pymc3 inference_library_version: 3.9.2 sampling_time: 18.096983432769775 - tuning_steps: 1000
  • created_at :
    2020-07-24T18:15:12.197697
    arviz_version :
    0.9.0
    inference_library :
    pymc3
    inference_library_version :
    3.9.2
    sampling_time :
    18.096983432769775
    tuning_steps :
    1000

  • - - + +
      @@ -6842,7 +6842,7 @@

      Example

      created_at: 2020-07-24T18:15:12.454586 arviz_version: 0.9.0 inference_library: pymc3 - inference_library_version: 3.9.2
    • g_coef
      PandasIndex
      PandasIndex(Index(['intercept', 'slope'], dtype='object', name='g_coef'))
  • created_at :
    2020-07-24T18:15:12.454586
    arviz_version :
    0.9.0
    inference_library :
    pymc3
    inference_library_version :
    3.9.2

  • - - + +
      @@ -7199,20 +7199,20 @@

      Example

      created_at: 2020-07-24T18:15:12.457652 arviz_version: 0.9.0 inference_library: pymc3 - inference_library_version: 3.9.2
  • created_at :
    2020-07-24T18:15:12.457652
    arviz_version :
    0.9.0
    inference_library :
    pymc3
    inference_library_version :
    3.9.2

  • - - + +
      @@ -7528,17 +7528,17 @@

      Example

      created_at: 2020-07-24T18:15:12.458415 arviz_version: 0.9.0 inference_library: pymc3 - inference_library_version: 3.9.2
  • created_at :
    2020-07-24T18:15:12.458415
    arviz_version :
    0.9.0
    inference_library :
    pymc3
    inference_library_version :
    3.9.2

  • - - + +
      @@ -7857,7 +7857,7 @@

      Example

      created_at: 2020-07-24T18:15:12.459832 arviz_version: 0.9.0 inference_library: pymc3 - inference_library_version: 3.9.2
  • created_at :
    2020-07-24T18:15:12.459832
    arviz_version :
    0.9.0
    inference_library :
    pymc3
    inference_library_version :
    3.9.2

    • @@ -8179,7 +8179,7 @@

    Example

    The function az.summary lets us look at some useful summary statistics, including \(\hat{R}\), divergent transitions and MCSE.

    The variable lp, which you can find in the group sample_stats is the model’s total log probability density. It’s not very meaningful on its own, but is useful for judging overall convergence. diverging counts the number of divergent transitions.

    -
    +
    az.summary(idata.sample_stats, var_names=["lp", "diverging"])
    /home/nicholas/courses/bayesian_statistics_for_computational_biology/.venv/lib/python3.12/site-packages/arviz/stats/diagnostics.py:592: RuntimeWarning: invalid value encountered in scalar divide
@@ -8238,7 +8238,7 @@ 
    Example

    Sometimes it’s useful to summarise individual parameters. This can be done by pointing az.summary at the group where the parameters of interest live. In this case the group is called posterior.

    -
    +
    az.summary(idata.posterior, var_names=["sigma", "g"])
    @@ -8305,7 +8305,7 @@

    Example

    The function az.loo performs approximate leave-one-out cross validation, which can be useful for evaluating how well the model might make predictions. Watch out for the warning column, which can tell you if the approximation is likely to be incorrect.

    -
    +
    az.loo(idata, var_name="y", pointwise=True)
    Computed from 2000 posterior samples and 919 observations log-likelihood matrix.
@@ -8324,7 +8324,7 @@ 
    Example

    The function az.compare is useful for comparing different out of sample log likelihood estimates.

    -
    +
    idata.log_likelihood["fake"] = xr.DataArray(
     # generate some fake log likelihoods
     np.random.normal(0, 2, [4, 500, 919]),
@@ -8372,7 +8372,7 @@ 
    Example
    
 -1027.176982
 26.817161
 0.000000
-0.984858
+0.984773
 28.848998
 0.000000
 False
@@ -8381,12 +8381,12 @@ 
    Example
    
 
 fake
 1
--1798.110835
-3627.506938
-770.933854
-0.015142
-3.461342
-28.940672
+-1797.875028
+3619.073150
+770.698047
+0.015227
+3.525172
+29.086009
 True
 log
 
diff --git a/materials/metropolis-hastings.qmd b/materials/metropolis-hastings.qmd
index e19c655..964ac66 100644
--- a/materials/metropolis-hastings.qmd
+++ b/materials/metropolis-hastings.qmd
@@ -235,6 +235,13 @@ Here is how you'd call the proposal function
 prop_mu, prop_sigma = proposal_multi(current_mu, current_sigma)
 ```
 
+the accept_reject probability should also be updated to account
+for the log-probability
+
+```
+accept_reject = np.exp(prob_prop - prob_current)
+```
+
 You 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.