Metropolis Hastings - Acceptance ratio, proposal and lkelihood - Cross Validated most recent 30 from stats.stackexchange.com 2019-07-16T17:18:01Z https://stats.stackexchange.com/feeds/question/378060 http://www.creativecommons.org/licenses/by-sa/3.0/rdf https://stats.stackexchange.com/q/378060 0 Metropolis Hastings - Acceptance ratio, proposal and lkelihood youpilat13 https://stats.stackexchange.com/users/226073 2018-11-21T04:04:28Z 2018-11-22T23:05:57Z <p>From a previous post :</p> <blockquote> <p>First to explain the MH algorithm, it's used to approximate numerically a <em>target distribution</em>, in this case <span class="math-container">$p(\theta|D)$</span>. At each stage of the algorithm:</p> <ol> <li>A value <span class="math-container">$\theta_{proposed}$</span> is proposed using the <em>jumping</em> or <em>proposal distribution</em>. </li> <li>An acceptance ratio is calculated, equal to <span class="math-container">$\frac{p(\theta_{proposed}|D)}{p(\theta_{current}|D)}$</span>. Because we cannot calculate <span class="math-container">$p(\theta|D)$</span> directly, we leverage the proportional expression of Bayes rule and calculate this quotient using the products of <em>likelihood</em> and <em>prior</em> corresponding to <span class="math-container">$\theta_{proposed}, \theta_{current}$</span>. That is, the acceptance ratio is:</li> </ol> <p><span class="math-container">$$\frac{p(D|\theta_{proposed})p(\theta_{proposed})}{p(D|\theta_{current})p(\theta_{current})}$$</span></p> <ol start="3"> <li>If this ratio exceeds one—intuitively, if the proposed parameter value is more likely given data and prior—we accept this proposal. If not, we accept it with probability equal to the ratio.</li> </ol> <p>As more and more values are sampled in this way, the trace of <span class="math-container">$\theta$</span> values more and more closely approximates the true distribution of <span class="math-container">$\theta$</span>.</p> </blockquote> <p><strong>Question 1)</strong> In my case, I have 6 parameters to estimate : Are <span class="math-container">$p(\theta_{\text{proposed}})$</span> and <span class="math-container">$p(\theta_{\text{current}})$</span> different ? especially if I am using uniform law.</p> <p><strong>Question 2)</strong> If yes, I could simply write :</p> <p><span class="math-container">$$\text{Acceptance}_{\text{Ratio}} = \dfrac{p(D|\theta_{\text{proposed}})}{p(D|\theta_{\text{current}})}$$</span></p> <p><strong>Question 3)</strong> If I use a normal distribution with variance and mean for each of 6 parametric, have I still the symetry of proposal ?</p> <p><strong>Question 4)</strong> How to express this acceptance ratio as a function of Likehood function or cost function used to find the 6 parameters ?</p> <p>My cost function is equal to :</p> <p><span class="math-container">$$J(\theta)=(y-H\,\theta)(y- H\,\theta)^T$$</span></p> <p>with <span class="math-container">$d = H\,\theta$</span> where <span class="math-container">$d$</span> is the response of system.</p> <p><strong>UPDATE 1:</strong> Here is below the Metropolis-Hastings I am using :</p> <pre><code> %%%%%% MCMC - Metropolis Algorithm %%%%%%% % Declare parameters of array nb = 1e5; % Metropolis loop % Index of main loop i = 1; % Reject number j = 0; % Accept number k = 0; % Starting values l1 = [10, 3]; l2 = [10, 3]; l3 = [10, 3]; l4 = [10, 3]; l5 = [10, 3]; l6 = [10, 3]; % Initial param array paramsArray = zeros(nb,6); paramsArray(1,1:6) = [l1(1) + l1(2)*rand(1), l2(1) + l2(2)*rand(1), ... l3(1) + l3(2)*rand(1), l4(1) + l4(2)*rand(1), ... l5(1) + l5(2)*rand(1), l6(1) + l6(2)*rand(1)]; % Init parameters vector pStart = paramsArray(1,1:6)'; % Standard deviation vectors for proposal distributon gamMatrix = diag([l1(2), l2(2), l3(2), l4(2), l5(2), l6(2)]); % Starting prior distribution f(p|d) w_x = Crit_J(pStart,D)*exp(-((pStart(1)-l1(1))^2/(2*l1(2)^2)+(pStart(2)-l2(1))^2/(2*l2(2)^2)+ ... (pStart(3)-l3(1))^2/(2*l3(2)^2)+(pStart(4)-l4(1))^2/(2*l4(2)^2)+ ... (pStart(5)-l5(1))^2/(2*l5(2)^2)+(pStart(6)-l6(1))^2/(2*l6(2)^2))); while (i &lt;= nb) if (i == 1) % First random ptest = pStart; else % Other random ptest = abs(pStart + gamMatrix*randn(6,1)); end % Upper term for acceptance ratio w_y = Crit_J(ptest,D)*exp(-((ptest(1)-l1(1))^2/(2*l1(2)^2)+(ptest(2)-l2(1))^2/(2*l2(2)^2)+ ... (ptest(3)-l3(1))^2/(2*l3(2)^2)+(ptest(4)-l4(1))^2/(2*l4(2)^2)+ ... (ptest(5)-l5(1))^2/(2*l5(2)^2)+(ptest(6)-l6(1))^2/(2*l6(2)^2))); % Generate u uniformly u = rand(1); % Ratio acceptance log_prob = log(w_y/w_x); % Test acceptation if (log(u) &lt; log_prob) % Assing new paramsArray paramsArray(i,1:6) = ptest(1:6); w_x = Crit_J(paramsArray(i,1:6),D)*exp(-((ptest(1)-l1(1))^2/(2*l1(2)^2)+(ptest(2)-l2(1))^2/(2*l2(2)^2)+ ... (ptest(3)-l3(1))^2/(2*l3(2)^2)+(ptest(4)-l4(1))^2/(2*l4(2)^2)+ ... (ptest(5)-l5(1))^2/(2*l5(2)^2)+(ptest(6)-l6(1))^2/(2*l6(2)^2))); i = i+1; k = k+1; else % Assing to previous if (i ~= 1) paramsArray(i,1:6) = paramsArray(i-1,1:6); i = i+1; j = j+1; end end end disp('acceptationt : ratio'); disp(k/nb) disp('reject : ratio'); disp(j/nb) % Display mean of different parameters disp('Parameters with Metropolis-Hastings :'); mean(paramsArray(:,1)) mean(paramsArray(:,2)) mean(paramsArray(:,3)) mean(paramsArray(:,4)) mean(paramsArray(:,5)) mean(paramsArray(:,6)) </code></pre> <p>and my cost function (assimilated to Likelihood function) :</p> <pre><code>% Function of cost function cost = Crit_J(p,D) % Compute the model corresponding to parameters p [R,C] = size(D); [Cols,Rows] = meshgrid(1:C,1:R); % Model Model = (1+((Rows-p(3)).^2+(Cols-p(4)).^2)/p(5)^2).^(-p(6)); model = Model(:); d = D(:); % Introduce H matrix H = [ model, ones(length(model),1)]; % Compute the cost function : taking absolute value cost = abs((d-H*[p(1),p(2)]')'*(d-H*[p(1),p(2)]')); end </code></pre> <p>If you could see the error ...</p> https://stats.stackexchange.com/questions/378060/-/378207#378207 3 Answer by Ben for Metropolis Hastings - Acceptance ratio, proposal and lkelihood Ben https://stats.stackexchange.com/users/173082 2018-11-22T02:21:52Z 2018-11-22T23:05:57Z <p>The problem with your description of the Metropolis-Hastings algorithm is that your notation does not distinguish between the probability densities in the actual problem you are trying to solve, and the <em>proposal density</em> used in the algorithm. Your notation also fails to capture the fact that we are trying to simulate from the posterior distribution, but we only have a kernel of this distribution. A better description of the algorithm, which makes these distinctions in the notation, is as follows:</p> <blockquote> <p>You start in a situation where you do not know the posterior <span class="math-container">$p(\theta|D)$</span>, but you do know a kernel of this distribution <span class="math-container">$K(\theta|D) \propto p(\theta|D)$</span>. You want to simulate values from the posterior. In the <strong>MH algorithm</strong> you start at an arbitrary parameter value <span class="math-container">$\theta_0$</span> and simulate using the following recursive scheme (which is a Markov chain):</p> <ol> <li>We generate a proposed value <span class="math-container">$\theta'_{t+1}$</span> from the <em>proposal density</em> <span class="math-container">$g(\theta'_{t+1}|\theta_{t})$</span>.</li> <li>For the proposed value, we define the acceptance ratio:</li> </ol> <p><span class="math-container">$$A(\theta'_{t+1} | \theta_t) \equiv \frac{K(\theta'_{t+1}|D)}{K(\theta_{t}|D)} \cdot \frac{g(\theta_{t} | \theta'_{t+1})}{g(\theta'_{t+1} | \theta_t)}.$$</span></p> <ol start="3"> <li>With probability <span class="math-container">$\min (A(\theta'_{t+1} | \theta_t), 1)$</span> we accept the proposed value and set <span class="math-container">$\theta_{t+1} = \theta'_{t+1}$</span>. Otherwise we reject the proposed value and set <span class="math-container">$\theta_{t+1} = \theta_{t}$</span>.</li> </ol> <p>It can be shown that this Markov chain has stationary distribution <span class="math-container">$p(\theta|D)$</span>, which is the posterior distribution of interest. Note that this is true even though the algorithm only uses a kernel of the distribution. We can therefore rely on the limiting properties of Markov chains to simulate from this posterior distribution. Usually this involves generating a small amount of 'burn-in' iterations followed by a series of auto-correlated simulations from the limiting stationary distribution. We can also rely on ergodic theorems to estimate the true posterior moments of functions of the parameter from the corresponding sample moments from the Markov chain.</p> <p><strong>Special case - symmetric proposal distribution:</strong> In many applications of the MH algorithm it is common to use a proposal density that is symmetric, in the sense that:</p> <p><span class="math-container">$$g(\theta'|\theta) = g(\theta|\theta') \quad \text{for all } \theta, \theta'.$$</span> (Note that a sufficient condition for this is that the density value depends on the parameters only through the norm <span class="math-container">$||\theta' - \theta||$</span>.) In this special case the acceptance ratio reduces to:</p> <p><span class="math-container">$$A(\theta'_{t+1} | \theta_t) \equiv \frac{K(\theta'_{t+1}|D)}{K(\theta_{t}|D)}.$$</span></p> </blockquote> <p>Now that we have a clearer explanation of the actual workings of the algorithm, I will try to answer your specific questions. (For consistency, I will translate your questions into notation that is consistent with my explanation of the algorithm.) Your Question 4 is unclear to me (there is no cost function in the algorithm so I don't know what you're referring to here), but I will answer the other three questions.</p> <hr> <blockquote> <p><strong>Question 1)</strong> Are <span class="math-container">$g(\theta'_{t+1} | \theta_t)$</span> and <span class="math-container">$g(\theta_t | \theta'_{t+1})$</span> different? What if I am using a uniform proposal distribution.</p> </blockquote> <p>In the case where you use a proposal distribution that is symmetric (in the sense described above) the two proposal densities (with the argument and conditioning parameter switched) will be the same. Symmetry occurs in the case where you use the uniform proposal density that is centred around the conditioning value:</p> <p><span class="math-container">$$g(\theta' | \theta) \propto \mathbb{I}(|\theta' - \theta| \leqslant t).$$</span></p> <p>In this case, switching the terms in the proposal density does not alter the value (i.e., they are not different). If you are using a uniform proposal density that is not centred around the conditioning value then this will not hold.</p> <blockquote> <p><strong>Question 2)</strong> If I am using a uniform proposal distribution, could I write the acceptance ratio without the ratio of proposal densities?</p> </blockquote> <p>Assuming your uniform proposal distribution is centred around the conditioning parameter (and thus symmetric in the above sense), yes you can.</p> <blockquote> <p><strong>Question 3)</strong> If I use a normal distribution, will I still have symmetry of the proposal distribution?</p> </blockquote> <p>Symmetry occurs if you use a normal distribution with mean equal to the conditioning parameter and variance independent of this parameter:</p> <p><span class="math-container">$$g(\theta' | \theta) = \text{N}(\theta' |\theta, \Sigma) \propto \exp \Big( -\frac{1}{2} (\theta' - \theta)^\text{T} \Sigma^{-1} (\theta' - \theta) \Big).$$</span></p>