Bayesian statistics: MCMC sampling technique
답변 (1개)
0 개 추천
Hi @Huthaifa,
To sample the posterior function using MCMC in MATLAB, you need to implement the Metropolis-Hastings algorithm. This algorithm is a popular choice for sampling from complex distributions where direct sampling is challenging. Below, I will provide a detailed explanation of the code, followed by the complete MATLAB implementation.
1. Define the Posterior Function: The posterior function is defined based on the likelihood and prior distributions. The likelihood is influenced by the number of positive and negative tests, while the prior reflects our initial beliefs about the parameters. 2. MCMC Setup: initialize parameters for the MCMC, including the number of samples, the proposal distribution, and the initial value. 3. Metropolis-Hastings Algorithm: This algorithm will be implemented to generate samples from the posterior distribution. It involves proposing a new sample and accepting or rejecting it based on the acceptance ratio. 4. Results Analysis: After generating the samples, compute the mean and standard deviation of the sampled values to summarize the posterior distribution.
%% Inputs % Likelihood mu_lampda_R = 1.2; % Average of lognormally distributed bias factors for the specific site. COV_lampda_R = 0.25; % Assumed COV (within site variability) for lampda R of the site. lampda_T = 1.0; % Minimum required bias factor represents a minimum acceptable testing load for the proof test. n = 3; % Total number of pile load tests conducted. m = 0; % Number of negative pile load tests. p = n - m; % Number of positive pile load tests.
% Prior mu_prime_lampda_R = 1.05; COV_Prime_lampda_R = 0.31;
%% Transfer from lognormal to normal space. % Likelihood sigma_ln_lampda_R = sqrt(log(1 + (COV_lampda_R)^2)); % Std dev of normally distributed resistance factor. mu_ln_lampda_R = log(mu_lampda_R) - 0.5 * (sigma_ln_lampda_R)^2; % Mean of normally distributed resistance factor for the site. ln_lampda_T = log(lampda_T);
% Prior sigma_prime_lampda_R = mu_prime_lampda_R * COV_Prime_lampda_R; % Std dev of lognormally distributed resistance factor. sigma_prime_ln_lampda_R = sqrt(log(1 + (COV_Prime_lampda_R)^2)); % Std dev of normally distributed resistance factor. mu_prime_ln_lampda_R = log(mu_prime_lampda_R) - 0.5 * (sigma_prime_ln_lampda_R)^2; % Mean of normally distributed resistance factor for the site. mu_prime_ln_mu = mu_prime_ln_lampda_R; % Mean for posterior. sigma_prime_ln_mu = sqrt((sigma_prime_ln_lampda_R)^2 - (sigma_ln_lampda_R)^2); % Std dev for posterior.
% Posterior function definition
posterior_fn = @(x) exp(-1 * (x - mu_prime_ln_mu)^2 / (2 *
sigma_prime_ln_mu^2)) * ...
(normcdf(-1 * ((ln_lampda_T - x) / sigma_ln_lampda_R)))^p * ...
((1 - normcdf(-1 * ((ln_lampda_T - x) / sigma_ln_lampda_R))))^m; %% MCMC Sampling num_samples = 10000; % Number of samples to draw samples = zeros(num_samples, 1); % Preallocate samples samples(1) = mu_prime_ln_mu; % Initial value acceptance_count = 0; % Count of accepted samples
% Proposal distribution parameters proposal_std = 0.1; % Standard deviation for the proposal distribution
for i = 2:num_samples % Propose a new sample proposed_sample = samples(i-1) + proposal_std * randn;
% Calculate acceptance ratio
acceptance_ratio = posterior_fn(proposed_sample) /
posterior_fn(samples(i-1)); % Accept or reject the proposed sample
if rand < acceptance_ratio
samples(i) = proposed_sample; % Accept the proposed sample
acceptance_count = acceptance_count + 1; % Increment acceptance
count
else
samples(i) = samples(i-1); % Reject and keep the previous sample
end
end%% Results mean_sample = mean(samples); % Mean of the samples std_sample = std(samples); % Standard deviation of the samples
% Display results
fprintf('Mean of the posterior samples: %.4f\n', mean_sample);
fprintf('Standard deviation of the posterior samples: %.4f\n', std_sample);
fprintf('Acceptance rate: %.2f%%\n', (acceptance_count / num_samples) *
100);
Please see attached.
The code begins by defining the necessary parameters for the likelihood and prior distributions. The posterior function is defined using the likelihood and prior, which is crucial for the MCMC sampling. The Metropolis-Hastings algorithm is implemented in a loop, where new samples are proposed and accepted based on the acceptance ratio. Finally, the mean and standard deviation of the sampled values are calculated and displayed, along with the acceptance rate of the MCMC process.
If you have any further questions or need additional modifications, feel free to ask!
댓글 수: 4
Hi @Huthaifa,
I have reviewed your comments. In MCMC methods, especially those using the Metropolis-Hastings algorithm, the choice of proposal distribution and its parameters is critical for achieving convergence to the target distribution (posterior). A low acceptance rate, such as 0.04%, suggests that proposed samples are often rejected, which can lead to poor exploration of the parameter space. When `COV_lampda_R` is greater than `COV_Prime_lampda_R`, it indicates that your likelihood is more variable than your prior beliefs, which can create difficulties in drawing meaningful samples. Here are some suggested solutions that can help.
1. Adjust Proposal Distribution: Increase the standard deviation of the proposal distribution (`proposal_std`). A larger standard deviation may help generate proposed samples that are further away from current samples, potentially increasing acceptance rates. Experiment with different distributions for proposals (e.g., Cauchy or uniform distributions) instead of a Gaussian distribution, as these can sometimes provide better exploration in high-dimensional spaces.
2. Tune Acceptance Ratio: Consider implementing adaptive MCMC techniques where you adjust `proposal_std` dynamically based on recent acceptance rates. For instance, if your acceptance rate falls below a certain threshold (like 20%), you might increase `proposal_std`.
3. Reassess Parameter Initialization: Ensure that your initial sample is reasonable and within a plausible range for your posterior distribution. Sometimes starting too far from where most of your posterior mass lies can lead to poor convergence.
4. Increase Sample Size: While you have set `num_samples = 10000`, increasing this number may provide more robust statistics, especially if you are seeing signs of non-convergence.
5. Examine Posterior Function: Double-check your implementation of the posterior function to ensure there are no logical errors or miscalculations, particularly in how likelihood and prior are combined.Implement diagnostic checks to visualize the posterior samples over iterations to identify if they are stuck in local minima.
6. Check for Infinite Values: Ensure that there are no divisions by zero or logarithms of zero in your calculations, as these can lead to undefined or infinite values.
Furthermore I have listed some additional insights that might help as well.
Convergence Diagnostics: Consider using convergence diagnostics like trace plots or Gelman-Rubin statistics to evaluate whether your chains have converged properly.
Alternative MCMC Methods: If issues persist, exploring other MCMC methods such as Hamiltonian Monte Carlo (HMC) or No-U-Turn Sampler (NUTS) could be beneficial as they often provide better convergence properties for complex distributions.
Parameter Sensitivity Analysis:Conduct sensitivity analyses to understand how changes in parameters like `COV_lampda_R` and `COV_Prime_lampda_R` affect convergence and sample quality.
The challenges you’re facing with MCMC sampling using Metropolis-Hastings may stem from an inappropriate proposal distribution or parameter settings that do not align well with your posterior landscape. By adjusting proposal parameters, reassessing initial values, and ensuring proper function definitions, you should see improvements in convergence and sampling efficiency.
If further issues arise or if specific adjustments yield unexpected results, please feel free to share additional details for more tailored guidance.
Hi @Huthaifa,
The Delayed Rejection Adaptive Metropolis (DRAM) algorithm is a sophisticated Markov Chain Monte Carlo (MCMC) method that enhances sampling efficiency. To implement DRAM in MATLAB, you can start with a basic structure that includes the proposal distribution, acceptance criteria, and the delayed rejection mechanism. Below is a simplified example to get you started:
function samples = dram(num_samples, initial_value) samples = zeros(num_samples, 1); samples(1) = initial_value; current_value = initial_value;
for i = 2:num_samples
% Propose a new value
proposed_value = current_value + randn;
% Simple Gaussian proposal % Calculate acceptance probability
acceptance_prob = min(1, exp(-0.5 * (proposed_value^2 -
current_value^2))); % Accept or reject the proposed value
if rand < acceptance_prob
current_value = proposed_value;
else
% Delayed rejection
proposed_value2 = current_value + randn; % Second proposal
acceptance_prob2 = min(1, exp(-0.5 * (proposed_value2^2 -
current_value^2)));
if rand < acceptance_prob2
current_value = proposed_value2;
end
end
samples(i) = current_value;
end
endThis code provides a basic framework for DRAM. You can expand upon it by incorporating more complex proposal distributions and tuning parameters to suit your specific application. As you gain more experience with MATLAB, you can refine and optimize this implementation further.
Hope this helps.
카테고리
도움말 센터 및 File Exchange에서 Stable Distribution에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
