montecarlo.m

蒙特卡洛方法

%this program is a monte carlo computer simulation that was
%used to generate figure 2.10a
%set seed of random number generator to initial value 
randn('seed',0);
%set up values of variance, data record length, and number 
%of realizations 
var=10;
N=10;
M=1000;
%dimension array of realizations
T=zeros(M,1);
%compute realization of the sample mean
for i=1:M
    x=sqrt(var)*randn(N,1);
    T(i)=mean(x);
end
%set number of values of gamma.
ngam=100;
%set up gamma array
gammamin=min(T);
gammamax=max(T);
gamdel=(gammamax-gammamin)/ngam;
gamma=[gammamin;gamdel;gammamax]';
%dimension P(the monte carlo estimate )and ptrue
%(the theoretical or true probability).
P=zeros(length(gamma),1);Ptrue=P;
%determine for each gamma how many realizations exceeded
%gamma(mgam) and use this to extimate the probability
for i=1:length(gamma)
    clear Mgam;
    Mgam=find(T>gamma(i));
    P(i)=length(Mgam)/M;
end
%compute the true probability;
% Ptrue=Q(gamma/(sqrt(var/N)));%计算高斯右尾概率
Ptrue=0.5*erfc((gamma/(sqrt(var/N)))/sqrt(2));
plot(gamma,P,'-',gamma,Ptrue,'--')
xlabel('gamma');
ylabel('P(T>gamma)');
grid