變化是生活的調味料

以下利用Matlab program來分別以Binomial RV及Bernoulli RV
%來試驗中央極限定理(Central Limit Theorem)
%Binomial RVs
%-------------MATLAB Script-------------------------------------
% file binomial_rv

n=input('Number of Bernoulli trials n=');
p=input('Parameter of Bernoulli RV p=');
m=input('Number of sample mean m=');

for k=1:m
i=0;
f=(1.0-p)^n;
q=f;
u=rand();
while f<=u
q=q*(p*(n-i)/((i+1.0)*(1.0-p)));
f=f+q;
i=i+1;
end
Binomial(k,1)=i;
end

mean=sum(Binomial)/m;
x=[1:m]';
Samples=[x Binomial];

% fprintf('\n\t Result of Binomial RVs\n')
% disp(Samples);

fprintf('The mean is = %f \n', mean);
---------------------------------------------------------------------------
>> binomial_rv
Number of Bernoulli trials n=6
Parameter of Bernoulli RV p=0.3
Number of sample mean m=10000
The mean is = 1.791700
>>
>> binomial_rv
Number of Bernoulli trials n=6
Parameter of Bernoulli RV p=0.3
Number of sample mean m=100000
The mean is = 1.801410
>>
The E[X]=np=1.8
As the number of samples approaches infinity, the sample mean approaches the E[X].
------------------------------------------------------

%以Bernoulli RV來測試CLT
%Bernoulli RVs
-------------MATLAB Script-------------------------------------
% file bernoulli_rv.m
n=input('Number of Bernoulli trials n=');
p=input('Parameter of Bernoulli RV p=');
B=rand(n,1)<=p;
mean=sum(B)/n;
x=[1:n]';
Samples=[x B];
% fprintf('\n\t Result of Binomial RVs\n')
% disp(Samples);
fprintf('The mean is = %f \n', mean);
----------------------------------------------------------------------
>> bernoulli_rv
Number of Bernoulli trials n=10000
Parameter of Bernoulli RV p=0.4
The mean is = 0.388700
>>
>> bernoulli_rv
Number of Bernoulli trials n=100000
Parameter of Bernoulli RV p=0.4
The mean is = 0.399340
>>
The E[X]=p=0.4.

It is very obvious that as the trials approaches infinity, the sample mean will approach E[x].

Write a program to generate Gaussian random variables. And then use this program to calculate E[X^2]. Using the moment generation function, we can get E[X^2]=μ2+σ2, and thus E[X^2]=0+1=1.
-------------MATLAB Script-------------------------------------
% file Gaussian_RV
u=input('mean of the gaussian RV u=');
var=input('variance of the gaussian RV var=');
m=input('Number of samples m=' );
temp=0;
for k=1:m
t=0;
if t==0;
r=2;
while r>1.0
v1=2.0*rand()-1.0;
v2=2.0*rand()-1.0;
r=v1*v1+v2*v2;
end
r=sqrt((-2.0*log(r))/r);
t=v2*r;
G(k,1)=u+v1*r*var;
else
temp=t
t=0.0;
G(k,1)=u+temp*var;
end
end
G2=G.*G;
mean=sum(G2)/m;
fprintf('The mean is = %f \n', mean);
------------------------------------------------------------------------
>> gaussian_rv
mean of the gaussian RV u=0
variance of the gaussian RV var=1
Number of samples m=10000
The mean is = 0.978683
>> gaussian_rv
mean of the gaussian RV u=0
variance of the gaussian RV var=1
Number of samples m=50000
The mean is = 0.997460
>>
-------------------------------------------------------------------------
We can see that the mean calculated using program is approaching the exact value of E[X^2] as the number of samples increases.

%(1)
%使用Matlab來產生Bernoulli Random Variable序列之和的樣本空間
%----------------------MATLAB Script-------------------------------------
% poisson_rvs
% Sequence of Poisson Random variables
n=input('Number of RVs n=');
m=input('Number of samples of the sum m=');
u=input('Parameter of Bernoulli RV u=');

% claculate the sum of n Poisson Rvs
for j=1:m
for k=1:n
i=0;
f=exp(-u);
p=f;
u=rand();
while f<=u
p=p*(u/(i+1.0));
f=f+p;
i=i+1;
end
P(k,j)=i;
end
end

for l=1:m
SUMOFPRV(l,1)=sum(P(:,l));
end

x=[1:m]';
Samples=[x SUMOFPRV];

fprintf('\n\tSome samples of the sum of these %d RVs',n)
fprintf('\n\t i \t S(i)\n');
disp(Samples);
------------------------------------------------------------------------
>> poisson_rvs
Number of RVs n=100
Number of samples of the sum m=10
Parameter of Poisson RV u=0.3

Some samples of the sum of these n RVs
i S(i)
1 63
2 46
3 55
4 37
5 60
6 51
7 47
8 48
9 47
10 34

%使用Matlab來產生Bernoulli Random Variable序列之和的樣本空間 
%------------------------------MATLAB Script-------------------------------------- 
% file bernoulli_rvs 
% Sequence of Bernoulli Random variables 
% X_i are Bernoulli RVs with i=1,...,n 
n=input('Number of RVs='); 
m=input('Number of samples of the sum=');
 p=input('Parameter of Bernoulli RV='); 

%Sn=X_1+X_2+...+X_100 
B=rand(n,m)<=p; 
for i=1:m 
SUMOFBRV(i,1)=sum(B(:,i)); 
end k=[1:m]'; 
Samples=[k SUMOFBRV]; 
fprintf('\n\tSome samples of the sum of these %d Bernoulli RVs',n);
 fprintf('\n\t i \t S(i)\n'); 
disp(Samples); 
-------------------------------------------------------------------------- 
>> bernoulli_rvs Number of RVs n=100 
Number of samples of the sum m=10
Parameter of Bernoulli RV p=0.3 
>>
 random_sum 
Some samples of the sum of these 100 RVs
 i S(i) 
1 25 
2 34 
3 36
4 36 
5 36 
6 32 
7 32 
8 29 
9 31 
10 39