Question

MONTE CARLO SIMULATION:

So I am trying to figure out a MATLAB code to generate a Monte Carlo simulation for my mean and standard deviation values. I have a material (fat) with mean and standard deviation values as well as mean and standard deviation values for its thickness. I need to figure out a code that will plot and generate a monte carlo simulation. Any help would be greatly appreciated. Thanks

Answer 1

% Imagine that we have obtained two samples with a correlation of r = 0.4

% for a sample size of 10. We do not know the population or have any other

% information.

%

% But we can assume that the true populations generating the two variables

% were normally distributed. We can then simulate a process in which we

% draw samples of size 10 from two independent Gaussian distributions and

% compute the correlation value observed in each sample.

%

% Such (Monte Carlo) simulations reveal that it is actually quite likely to

% obtain a correlation of r = 0.4 for a sample size of 10, even when the

% underlying distributions are independent. Thus, we should not believe the

% claim.

% Empirical (measured) correlation.

r_empirical = .4;

% We simulate two independent and normally distributed variables.

N = 10000; % We set the number of montecarlo simulations

n = 10; % The sample size

s1_sim = randn(N,n); % siulated distribution of samples

s2_sim = randn(N,n); % simulated distribution of samples

% compute the correlation coefficient

for ii = 1:N

  % Note, here we are using matlab's corr.m to compute the pearon

  % correlation instead of computing it manually.

  r_dist(ii) = corr(s1_sim(ii,:)',s2_sim(ii,:)','type' , 'Pearson' );

end

% We plot the distribution of correlation coefficients obtained this way.

figure('name','Test Correlation (Monte Carlo Method)','color','w');

% Let's compute the histogram of the distribution and trasform the

% distribution into a probabilities.

%

% We choose the sampling spacing between the bins.

binWidth = .01;

% We make bins that cover the full range fo values in P and are sampled at

% the resolution we decided using binWidth.

bins = min(r_dist):binWidth:max(r_dist);

% We now make a histogram of the distribution of these differences.

% We scale the dstributions properly to transform into probabilites.

[nn, xx] = hist(r_dist,bins);

bar(xx,nn / (N * binWidth),'k');

title('Distribution of the Montecarlo-simulated correlation coefficients','FontSize',14)

ylabel('Probability of occurrence','FontSize',14)

xlabel('Corr(s1,s2)','FontSize',14)

hold on

set(gca,'tickdir','out', 'box','off','FontSize',14)

% Get the limits of the y axis from the plot. We will use these for the

% next plots.

y = get(gca,'yLim');

% We can now show the empirical difference we simulated (d_empirical). The

% height of the distribution at the value of d_empirical is the probability

% that that difference comes from samples drawn from the same distribution.

plot([r_empirical,r_empirical],y*.99,'r-','lineWidth',2)

% We can obtain the probability of H0 being false by counting the number of

% randomly obtained correlation coefficient values that are larger than the

% empirical value and divide this by the total number of simulations that

% were run. Because the value can be larger or smaller than the mean we

% coupute the absolute value of the ditribution and r values.

p = sum(abs(r_empirical) < abs(r_dist))/N;

text(-0.8,y(2)*.95,sprintf('H0 is false with %2.4f probability.',p),'FontSize',14)