matlab 按照概率生成数字_matlab生成的随机数是真正随机的吗？

. `9 Y) n7 _7 N5 G1 K( Q

随机数序列在数值分析和概率统计中占有非常重要的地位，因为使用蒙特卡罗模拟方法的前提就是要求很多足够多的，真正的随机数。matlab是基于某种算法，通过rand函数来产生随机数的。从随机数的定义看，rand函数产生的序列不是随机数，是伪随机数。但我们在使用蒙特卡罗模拟方法算法时，不可能成千上万次的去投掷硬币来产生随机数，所以要考察matlab产生的伪随机数能不能当随机数使用。6 t9 b( t" E: F/ b3 l

考察的方法是：2 v9 T9 i" [9 M

1：利用rand函数产生200个伪随机数，分别统计出有多少个伪随机数小数点后的5位数字中奇数出现的频数为0，有多少个伪随机数小数点后的5位数字中奇数出现的频数为1，有多少个伪随机数小数点后的5位数字中奇数出现的频数为2，一直统计到有多少个伪随机数小数点后的5位数字中奇数出现的频数为5。 这样就获得一个包含6个元素的行向量s，其元素依次为小数点后5位包含奇数个数为0或1或2……或5的伪随机数的个数。) K) } J* c. v& E8 t) Y$ E

2：给出两组源于实际观测的数据。一组是记录某个医院相继出生的1000个婴儿的随机序列，记5个连续出生的婴儿为一组，这样共有200组，统计分别含有0，1，2……5个男婴的组数，获得向量a；另一组是从一个装有500黑球和500个白球的口袋里，每次有放回的抽取一个球，供抽取1000次，记连续抽样5次为一组，这样共有200组，统计分别含有0，1，2……5个白球的组数，获得向量b。; l4 W% A: [$ C. u# _9 S/ C. w/ Q

3：对向量a，b，s进行自由度为5卡方检验，分别获得以向量a，b，s为代表的三组数据的χ2值。& y5 O/ G6 @% M$ f7 u" _5 s

4：将上述步骤重复1000次，每次向量a，b都是不变的，但每次的s向量都不同。

7 P$ R% l& E' - e

5：计算1000组s对应的1000个χ2值的最大值，最小值，平均值，对1000组a，b同样如此。: j4 I. V/ A" `( e8 P& I4 W( o

6：如果假设显著水平为0.05，那么自由度为5的χ2分布临界值是11.1，所以还要计算1000组s对应的1000个χ2值中大于11.1的数值所占的分率。/ z% R7 g+ |( N

总结果如下：

$ ? R: E c h( o- U. L1 r, S

a b s(基于matlab)% ]! x* f# l5 t$ z

平均值：2.2240 5.0400 5.0038

( v3 K' x* ?9 e: i( y4 A$ O

最大值：2.2240 5.0400 19.3760

. d. C5 Q$ c* v& V9 B

最小值： 2.2240 5.0400 0.2560

% h# {$ c: X4 u! R" x/ Q

1000个χ2值中大于11.1的数值所占的分率7 v3 @! G1 d& r9 R' Z( j8 p

x = 0.0520

5 l) a; h: q4 ~' p' e. H T

从平均值的计算结果看，matlab产生的伪随机数的随机程度和从口袋摸球相当，所以随机性满足要求。* I" u. i4 ~- S

从最大值结果看，基于matlab的伪随机序列产生的χ2值最大达到19.3760，大于显著水平0.05，自由度为5的χ2分布的临界值，似乎有些序列不够随机。但考虑到χ2分布中，总有0.05的概率，使得χ2值大于11.1，所以验算了基于matlab的伪随机序列产生的χ2值中大于11.1的数值占的分率，这个分率是0.0520，非常接近0.05。

$ `) A1 n) b/ }& D8 m& R

所以，matlab产生的伪随机序列可以作为真正的随机序列使用。8 F8 k0 O1 R6 [; Q4 T; O

matlab程序如下：

7 U& 4 {; f% T% L8 p- g" z9 J9 y e6 }% ^5 f H# I Y

clear

rr=[];

for l=1:1000

p=[];

rand('seed',prod(clock))

r=fix(rand(200)*100000);

for i=1:length(r)

m=r(i); s=[];

for j=1:5

s=[s rem(m,10)];

m=round((m-rem(m,10))/10);

end

p=[p;s];

end

s=zeros(1,6);

for i=1:size(p,1)

k=length(find(rem(p(i,

,2*ones(1,size(p,2)))));

s(k+1)=s(k+1)+1;

end

a=[5 27 64 65 30 9];

b=[4 34 65 70 22 5];

c=[];

for i=0:5

p=combine_m(5,i)*0.5^i*0.5^(5-i);

c=[c 200*p];

end

ka2=sum([([a;b;s]-[c;c;c]).^2./[c;c;c]]');

rr=[rr;ka2];

end

ave=mean(rr)

mx=max(rr)

mi=min(rr)

x=length(find(rr(:,3)>11.1))/length(rr(:,3))