The textbook example of the Inverse Transform Method is the exponential distribution. The set of pseudorandom numbers given by<\/p>\n\n\n\n
\\[ x = -\\ln (1-y) \\; , \\]<\/p>\n\n\n\n
where $y$ is uniformly distributed on the interval $[0,1]$, is guaranteed to be distributed exponentially. <\/p>\n\n\n\n
If no such relationship can be found, the Inverse Transform Method fails to be useful. <\/p>\n\n\n\n
The Metropolis-Hastings method is capable to overcoming this hurdle at the expense of being computationally much more expensive. The algorithm works as follows:<\/p>\n\n\n\n
- Setup
- Select the distribution $f(x)$ from which a set random numbers are desired<\/li><\/ul>
- Arbitrarily assign a value to a random realization $x_{start}$ of the desired distribution consistent with its range<\/li><\/ul>
- Select a simple distribution $g(x|x_{cur})$, centered on the current value of the realization from which to draw candidate changes<\/li><\/ul><\/li>
- Execution
- Draw a change from $g(x|x_{cur})$<\/li><\/ul>
- Create a candidate new realization $x_{cand}$<\/li><\/ul>
- Form the ratio $R = f(x_{cand})\/f(x_{cur})$<\/li><\/ul>
- Make a decision as the keeping or rejecting $x_{cand}$ as the new state to transition to according to the Metropolis-Hastings rule
- If $R \\geq 1$ then keep the change and let $x_{cur} = x_{cand}$<\/li><\/ul>
- If $R < 1$ then draw a random number $q$ from the unit uniform distribution and make the transition if $q < R$<\/li><\/ul><\/li><\/ul><\/li><\/ul>\n\n\n\n
The bulk of this algorithm is very straightforward and simple. The only confusing part is the Metropolis-Hastings decision rule in the last step and why it requires the various pieces above to work. The next post will be dedicated to a detailed discussion of that rule. For now, we\u2019ll just show the algorithm in action with various forms of $f(x)$ which are otherwise intractable.<\/p>\n\n\n\n
The python code for implementing this algorithm is relatively simple. Below is a snippet from the implementation that produced the numerical results that follow. The key points to note are:<\/p>\n\n\n\n
- the numpy unit uniform distribution ($U(0,1)$) random<\/em><\/strong> was used for both $g(x|x_{cur})$ and the Metropolis-Hastings decision rule<\/li>
- N_trials<\/em><\/strong> is the number of times the decision loop was executed but is not the number of resulting accepted or kept random numbers, which is necessarily a smaller number<\/li>
- delta<\/em><\/strong> is the parameter that governs the size of the candidate change (i.e., $g(x|x_{cur}) = U(x_{cur}-\\delta,x_{cur}+\\delta)$<\/li>
- n_burn<\/em><\/strong> is the number of initial accepted changes to ignore or throw away (more on this next post)<\/li>
- x_0 <\/em><\/strong>is the arbitrarily assigned starting realization value<\/li>
- prob<\/em><\/strong> is a pointer to the function $f(x)$<\/li><\/ul>\n\n\n\n
In addition to these input parameters, the value $f = N_{kept}\/N_{trials}$ was reported. The rule of thumb for Metropolis-Hastings is to have a value $0.25 \\leq f \\leq 0.33$.<\/p>\n\n\n\n
\nx_lim = 2.8\nN_trials = 1_000_000\ndelta = 7.0\nDs = (2*random(N_trials) - 1.0)*delta\nRs = random(N_trials)\nn_accept = 0\nn_burn = 1_000\nx_0 = -2.0\nx_cur = x_0\nprob = pos_tanh\np_cur = prob(x_cur,-x_lim,x_lim)\n \nx = np.arange(-3,3,0.01)\nz = np.array([prob(xval,-x_lim,x_lim) for xval in x])\n \nkept = []\nfor d,r in zip(Ds,Rs):\n x_trial = x_cur + d\n p_trial = prob(x_trial,-x_lim,x_lim)\n w = p_trial\/p_cur\n if w >= 1 or w > r:\n x_cur = x_trial\n n_accept += 1\n if n_accept > n_burn:\n kept.append(x_cur)<\/pre><\/div>\n\n\n\n
Two arbitrary and difficult functions were used for testing the algorithm: $f_1(x) = |tanh(x)|\/N_1$ and $f_2(x) = \\cos^6 (x)\/N_2$. $N_1$ and $N_2$ are normalization constants obtained by numerically integrating over the range $[-2.8,2.8]$, which was chosen for strictly aesthetic reasons (i.e., the resulting plots looked nice).<\/p>\n\n\n\n
The two parameters that had the most impact of the performance of the algorithm were the values of delta<\/em><\/strong> and x_0<\/em><\/strong>. The metric chosen for the \u2018goodness\u2019 of the algorithm is simply a visual inspection of how well the resulting distributions matched the normalized probability density functions $f_1(x)$ and $f_2(x)$.<\/p>\n\n\n\n
Case 1 – $f_1(x) = |tanh(x)|\/N_1$<\/h2>\n\n\n\n
This distribution was chosen since it has a sharp cusp in the center and two fairly flat and broad wings. The best results came from a starting point $x_0 = -2.3$ over in the left wing (the right wing also worked equally well) and a fairly large value of $\\delta=7.0$ was chosen such that $f \\approx 1\/3$. The agreement is quite good.<\/p>\n\n\n\n
![\"\"](\"http:\/\/aristotle2digital.blogwyrm.com\/wp-content\/uploads\/2022\/05\/pos_tanh_x0_left_delta_big.png\")
<\/figure><\/div>\n\n\n\nA smaller value of $\\delta$ lead to more samples being collected ($f=0.85$) but the realizations couldn\u2019t \u2018shake off the memory\u2019 that they started over in the left wing and the results were not very good, although the two-wing structure can still be seen.<\/p>\n\n\n\n
![\"\"](\"http:\/\/aristotle2digital.blogwyrm.com\/wp-content\/uploads\/2022\/05\/pos_tanh_x0_left_delta_small.png\")
<\/figure><\/div>\n\n\n\nStarting off with $x_0 = 0$ resulted in a realization of random numbers that looks far more like the uniform distribution than the distribution we were hunting.<\/p>\n\n\n\n
![\"\"](\"http:\/\/aristotle2digital.blogwyrm.com\/wp-content\/uploads\/2022\/05\/pos_tanh_x0_middle_delta_big.png\")
<\/figure>\n\n\n\nThis result persists with smaller values of $\\delta$.<\/p>\n\n\n\n
Case 2 – $f_2(x) = \\cos^6(x)\/N_2$<\/h2>\n\n\n\n
In this case, the distribution has a large central hump or mode and two truncated wings. The best results came from a starting value of $x_0$ at or near zero and a large value of $\\delta$.<\/p>\n\n\n\n
![\"\"](\"http:\/\/aristotle2digital.blogwyrm.com\/wp-content\/uploads\/2022\/05\/cos_pos_x0_middle_delta_big.png\")
<\/figure><\/div>\n\n\n\nLowering the value of $\\delta$ leads to the algorithm being trapped in the central peak with no way to cross into the wings.<\/p>\n\n\n\n
![\"\"](\"http:\/\/aristotle2digital.blogwyrm.com\/wp-content\/uploads\/2022\/05\/cos_pos_x0_middle_delta_small.png\")
<\/figure><\/div>\n\n\n\nFinally, starting in the wings with a large value of $\\delta$ leads to a realization that looks like it has the general idea as to the shape of the distribution but without the ability to capture the fine detail.<\/p>\n\n\n\n
![\"\"](\"http:\/\/aristotle2digital.blogwyrm.com\/wp-content\/uploads\/2022\/05\/cos_pos_x0_right_delta_big.png\")
<\/figure><\/div>\n\n\n\nThis behavior persists even when the number of trials is increased by a factor of 10 and no immediately satisfactory explanation exists.<\/p>\n","protected":false},"excerpt":{"rendered":"
The last two columns explored the notion for transition probabilities, which are conditional probabilities that express the probability of a new state being realized some system given the current state… Read more ><\/a><\/p>\n","protected":false},"author":2,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[1],"tags":[],"class_list":["post-1467","post","type-post","status-publish","format-standard","hentry","category-uncategorized"],"_links":{"self":[{"href":"https:\/\/aristotle2digital.blogwyrm.com\/index.php?rest_route=\/wp\/v2\/posts\/1467","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/aristotle2digital.blogwyrm.com\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/aristotle2digital.blogwyrm.com\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/aristotle2digital.blogwyrm.com\/index.php?rest_route=\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/aristotle2digital.blogwyrm.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=1467"}],"version-history":[{"count":1,"href":"https:\/\/aristotle2digital.blogwyrm.com\/index.php?rest_route=\/wp\/v2\/posts\/1467\/revisions"}],"predecessor-version":[{"id":1481,"href":"https:\/\/aristotle2digital.blogwyrm.com\/index.php?rest_route=\/wp\/v2\/posts\/1467\/revisions\/1481"}],"wp:attachment":[{"href":"https:\/\/aristotle2digital.blogwyrm.com\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=1467"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/aristotle2digital.blogwyrm.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=1467"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/aristotle2digital.blogwyrm.com\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=1467"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}
- N_trials<\/em><\/strong> is the number of times the decision loop was executed but is not the number of resulting accepted or kept random numbers, which is necessarily a smaller number<\/li>
- the numpy unit uniform distribution ($U(0,1)$) random<\/em><\/strong> was used for both $g(x|x_{cur})$ and the Metropolis-Hastings decision rule<\/li>
- If $R < 1$ then draw a random number $q$ from the unit uniform distribution and make the transition if $q < R$<\/li><\/ul><\/li><\/ul><\/li><\/ul>\n\n\n\n
- If $R \\geq 1$ then keep the change and let $x_{cur} = x_{cand}$<\/li><\/ul>
- Make a decision as the keeping or rejecting $x_{cand}$ as the new state to transition to according to the Metropolis-Hastings rule
- Form the ratio $R = f(x_{cand})\/f(x_{cur})$<\/li><\/ul>
- Create a candidate new realization $x_{cand}$<\/li><\/ul>
- Draw a change from $g(x|x_{cur})$<\/li><\/ul>
- Arbitrarily assign a value to a random realization $x_{start}$ of the desired distribution consistent with its range<\/li><\/ul>
- Select the distribution $f(x)$ from which a set random numbers are desired<\/li><\/ul>