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!
