第8章 概率论与数理统计问题的求解 概率分布与伪随机数生成 统计量分析 数理统计分析方法及计算机实现 统计假设检验 方差分析及计算机求解.

Presentation on theme: "第8章 概率论与数理统计问题的求解 概率分布与伪随机数生成 统计量分析 数理统计分析方法及计算机实现 统计假设检验 方差分析及计算机求解."— Presentation transcript:

1 第8章 概率论与数理统计问题的求解 概率分布与伪随机数生成 统计量分析 数理统计分析方法及计算机实现 统计假设检验 方差分析及计算机求解

2 8.1概率分布与伪随机数生成 8.1.1 概率密度函数与分布函数概述
8.1概率分布与伪随机数生成 概率密度函数与分布函数概述

3 通用函数计算概率密度函数值 函数 pdf 格式 P=pdf(‘name’，K，A) P=pdf(‘name’，K，A，B) P=pdf(‘name’，K，A，B，C) 说明 返回在X=K处、参数为A、B、C的概率密度值，对于不同的分布，参数个数是不同；name为分布函数名。 例如二项分布：设一次试验，事件Y发生的概率为p，那么，在n次独立重复试验中，事件Y恰好发生K次的概率P_K为：P_K=P{X=K}=pdf('bino'，K，n，p)

4 例: 计算正态分布N（0，1）的随机变量X在点0.6578的密度函数值。
解： >> pdf('norm',0.6578,0,1) ans = 0.3213 例:自由度为8的卡方分布，在点2.18处的密度函数值。 解： >> pdf('chi2',2.18,8) 0.0363

5 随机变量的累积概率值(分布函数值) 通用函数cdf用来计算随机变量的概率之和（累积概率值） 函数 cdf 格式 cdf(‘name’，K，A) cdf(‘name’，K，A，B) cdf(‘name’，K，A，B，C) 说明 返回以name为分布、随机变量X≤K的概率之和的累积概率值，name为分布函数名.

6 例: 求标准正态分布随机变量X落在区间(-∞，0.4)内的概率。
解：>> cdf('norm',0.4,0,1) ans = 0.6554 例：求自由度为16的卡方分布随机变量落在[0，6.91]内的概率。 解：>> cdf('chi2',6.91,16) 0.0250

7 随机变量的逆累积分布函数 MATLAB中的逆累积分布函数是已知，求x。 命令 icdf 计算逆累积分布函数 格式 icdf(‘name’，F，A) icdf(‘name’，F，A，B) icdf(‘name’，F，A，B，C) 说明 返回分布为name，参数为a1,a2,a3，累积概率值为P的临界值，这里name与前面相同。 如果F= cdf(‘name’，X，A，B，C) ， 则 X = icdf(‘name’，F，A，B，C)

8 例:在标准正态分布表中，若已知F=0.6554,求X 解： >> icdf('norm',0.6554,0,1) ans = 0.3999 例：公共汽车门的高度是按成年男子与车门顶碰头的机会不超过1%设计的。设男子身高X（单位：cm）服从正态分布N（175，6），求车门的最低高度。 解：设h为车门高度，X为身高。 求满足条件 F{X>h}<=0.99，即 F{X<h}>=0.01故 >> h=icdf('norm',0.99, 175, 6) h =

9 8.1.2 常见分布的概率密度函数与分布函数 8.1.2.1 Poisson分布
其要求x是正整数。

10 其中：x为选定的一组横坐标向量， y为x各点处的概率密度函数值。

11 例：绘制 l =1,2,5,10 时 Poisson 分布的概率密度函数与概率分布函数曲线。
>> x=[0:15]'; y1=[]; y2=[]; lam1=[1,2,5,10]; >> for i=1:length(lam1) y1=[y1,poisspdf(x,lam1(i))]; y2=[y2,poisscdf(x,lam1(i))]; end >> plot(x,y1), figure; plot(x,y2)

12 正态分布 正态分布的概率密度函数为:

13 例： >> x=[-5:.02:5]'; y1=[]; y2=[];
>> mu1=[-1,0,0,0,1]; sig1=[1,0.1,1,10,1]; sig1=sqrt(sig1); >> for i=1:length(mu1) y1=[y1,normpdf(x,mu1(i),sig1(i))]; y2=[y2,normcdf(x,mu1(i),sig1(i))]; end >> plot(x,y1), figure; plot(x,y2)

14 分布

15 例：绘制 为（1,1),(1,0.5),(2,1),(1,2),(3,1)时 >> x=[-0.5:.02:5]‘; ％x=[-eps:-0.02:-0.5,0:0.02:5]; x=sort(x’);替代 >> y1=[]; y2=[]; a1=[1,1,2,1,3]; lam1=[1,0.5,1,2,1]; >> for i=1:length(a1) y1=[y1,gampdf(x,a1(i),lam1(i))]; y2=[y2,gamcdf(x,a1(i),lam1(i))]; end >> plot(x,y1), figure; plot(x,y2)

16 分布（卡方分布） 其为一特殊的 分布 ，a=k/2, l =1/2。

17 例： >> x=[-eps:-0.02:-0.5,0:0.02:2]; x=sort(x');
>> k1=[1,2,3,4,5]; y1=[]; y2=[]; >> for i=1:length(k1) y1=[y1,chi2pdf(x,k1(i))]; y2=[y2,chi2cdf(x,k1(i))];end >> plot(x,y1), figure; plot(x,y2)

18 分布 概率密度函数为: 其为参数k的函数，且k为正整数。

19 例： >> x=[-5:0.02:5]'; k1=[1,2,5,10]; y1=[]; y2=[];
>> for i=1:length(k1) y1=[y1,tpdf(x,k1(i))]; y2=[y2,tcdf(x,k1(i))]; end >> plot(x,y1), figure; plot(x,y2)

20 Rayleigh分布

21 例： >> x=[-eps:-0.02:-0.5,0:0.02:5]; x=sort(x');
>> b1=[.5,1,3,5]; y1=[]; y2=[]; >> for i=1:length(b1) y1=[y1,raylpdf(x,b1(i))]; y2=[y2,raylcdf(x,b1(i))]; end >> plot(x,y1), figure; plot(x,y2)

22 F 分布 其为参数p,q的函数，且p,q均为正整数。

23 例：分别绘制（p,q）为（1,1),(2,1),(3,1)(3,2),(4,1)时F分布的概率密度函数与分布函数曲线。
>> x=[-eps:-0.02:-0.5,0:0.02:1]; x=sort(x'); >> p1=[ ]; q1=[ ]; y1=[]; y2=[]; >> for i=1:length(p1) y1=[y1,fpdf(x,p1(i),q1(i))]; y2=[y2,fcdf(x,p1(i),q1(i))]; end >> plot(x,y1), figure; plot(x,y2)

24 8.1.3 概率问题的求解 图4-9

25 例： >> b=1; p1=raylcdf(0.2,b); p2=raylcdf(2,b); P1=p2-p1 P1 = 0.8449 >> p1=raylcdf(1,b); P2=1-p1 P2 = 0.6065

26 例： >> syms x y; f=x^2+x*y/3; >> P=int(int(f,x,0,1/2),y,0,1/2) P = 5/192 >> syms x y; f=x^2+x*y/3; P=int(int(f,x,0,1),y,0,2) 1

27 8.1.4 随机数与伪随机数

28

29 例： >> b=1; p=raylrnd(1,30000,1);
>> xx=0:.1:4; yy=hist(p,xx); ％ hist()找出随机数落入各个子区间的点个数，并由之拟合出生成数据的概率密度。 >>yy=yy/(30000*0.1); >> bar(xx,yy), >> y=raylpdf(xx,1); >> line(xx,y)

30 8.2 统计量分析 随机变量的均值与方差

31 例： 均值 >> syms x; syms a lam positive
>> p=lam^a*x^(a-1)/gamma(a)*exp(-lam*x); >> m=int(x*p,x,0,inf) m = 1/lam*a 方差 >> s=simple(int((x-1/lam*a)^2*p,x,0,inf)) s = a/lam^2

32 已知一组随机变量样本数据构成的向量: 求该向量各个元素的均值、方差和标准差、中位数median

33 例：生成一组 30000 个正态分布随机数，使其均值为 0. 5，标准差为1
例：生成一组 个正态分布随机数，使其均值为 0.5，标准差为1.5，分析数据实际的均值、方差和标准差，如果减小随机变量个数，会有什么结果？ >> p=normrnd(0.5,1.5,30000,1);[mean(p),var(p),std(p)] ans = 300个随机数 >> p=normrnd(0.5,1.5,300,1);[mean(p),var(p),std(p)] ％可见在进行较精确的统计分析时不能选择太小的样本点。

34 例： >> [m,s]=raylstat(0.45) m = 0.5640 s = 0.0869

35 8.2.2 随机变量的矩

36 例： 求解原点矩 >> syms x; syms a lam positive; p=lam^a*x^(a-1)/gamma(a)*exp(-lam*x); >> for n=1:5, m=int(x^n*p,x,0,inf), end m = 1/lam*a 1/lam^2*a*(a+1) 1/lam^3*a*(a+1)*(a+2) 1/lam^4*a*(a+1)*(a+2)*(a+3) 1/lam^5*a*(a+1)*(a+2)*(a+3)*(a+4) ％有规律

37 >> syms n; m=simple(int((x)^n*p,x,0,inf)) ％直接求出
lam^(-n)*gamma(n+a)/gamma(a) >> for n=1:6, s=simple(int((x-1/lam*a)^n*p,x,0,inf)), end ％中心距 s = a/lam^2 2*a/lam^3 3*a*(a+2)/lam^4 4*a*(5*a+6)/lam^5 5*a*(3*a^2+26*a+24)/lam^6 ％好像无规律

38

39 例：考虑前面的随机数，可以用下面的语句得出随机数的各阶矩。
>> A=[]; B=[]; p=normrnd(0.5,1.5,30000,1); n=1:5; >> for r=n, A=[A, sum(p.^r)/length(p)]; B=[B,moment(p,r)]; end >> A,B A = B =

40 求各阶距的理论值： >> syms x; A1=[]; B1=[]; p=1/(sqrt(2*pi)*1.5)*exp(-(x-0.5)^2/(2*1.5^2)); >> for i=1:5 A1=[A1,vpa(int(x^i*p,x,-inf,inf),12)]; B1=[B1,vpa(int((x-0.5)^i*p,x,-inf,inf),12)]; end >> A1, B1 A1 = [ , , , , ] B1 = [ , , , , ]

41 8.2.3 多变量随机数的协方差分析

42

43 例： >> p=randn(30000,4); cov(p) ans =

44 8.2.4 多变量正态分布的联合概率 密度即分布函数

45 例： >> mu1=[-1,2]; Sigma2=[1 1; 1 3]; % 输入均值向量和协方差矩阵
>> [X,Y]=meshgrid(-3:0.1:1,-2:0.1:4); xy=[X(:) Y(:)]; % 产生网格数据并处理(两列2501*2 ） >> p=mvnpdf(xy,mu1,Sigma2); % 求取联合概率密度 >> P=reshape(p,size(X)); ％Change size(2501*1—61*41) >> surf(X,Y,P)

46 对协方差矩阵进行处理，可计算出新的联合概率密度函数。
>> Sigma2=diag(diag(Sigma2)); % 消除协方差矩阵的非对角元素 >> p=mvnpdf(xy,mu1,Sigma2); P=reshape(p,size(X)); surf(X,Y,P) R为m行n列。

47 例： >> mu1=[-1,2]; Sigma2=[1 1; 1 3];
>> R1=mvnrnd(mu1,Sigma2,2000); plot(R1(:,1),R1(:,2),'o') >> Sigma2=diag(diag(Sigma2)); figure; >> R2=mvnrnd(mu1,Sigma2,2000); plot(R2(:,1),R2(:,2),'o')

48 8.3数理统计分析方法及计算机实现 8.3.1 参数估计与区间估计
8.3数理统计分析方法及计算机实现 参数估计与区间估计 无论总体X的分布函数 的类型已知或未知，我们总是需要去估计某些未知参数或数字特征，这就是参数估计问题.即参数估计就是从样本 出发，构造一些统计量 （i=1，2，…，k）去估计总体X中的某些参数（或数字特征） （i=1，2，…，k）.这样的统计量称为估计量.

49 1、点估计：构造 的函数 作为参数 的点估计量， 称统计量 为总体X参数 的点估计量. 2.   区间估计：构造两个函数 (X1，X2，…， Xn）和 (X1，X2，…， Xn）做成区间，把 这 作为参数 的区间估计.

50 区间估计的求法 设总体X的分布中含有未知参数 ，若对于给定的概率 ，存在两个统计量 (X1，X2，…，Xn）和 (X1，X2，…，Xn）,使得 则称随机区间 为参数 的置信水平为 的置信区间,称 为置信下限,称 为置信上限.

51 由极大拟然法估计出该分布的均值、方差 及其置信区间。置信度越大，得出的置信区间越小，即得出的结果越接近于真值。
还有gamfit(), raylfit(), poissfit() ，unifit()（均匀分布） 等参数估计函数

52 例： >> p=gamrnd(1.5,3,30000,1); Pv=[0.9,0.92,0.95,0.98]; A=[]; >> for i=1:length(Pv) [a,b]=gamfit(p,Pv(i)); A=[A; Pv(i),a(1),b(:,1)',a(2),b(:,2)'] end >> A A =

53 >> num=[300,3000,30000,300000,3000000]; A=[];
>> for i=1:length(num) p=gamrnd(1.5,3,num(i),1); [a,b]=gamfit(p,0.95); A=[A;num(i),a(1),b(:,1)',a(2),b(:,2)']; end >> A(:,[2,3,4,5,6,7]) ans = 要达到参数估计效果良好，随机数不能选得太少，也不能选得太多，此例中为30000为好。

54 8.3.2 多元线性回归与区间估计

55

56

57 例： >> a=[1 -1.232 2.23 2 4 3.792]'; X=randn(120,6); y=X*a;
>> a1=inv(X'*X)*X'*y;a1' ans = >> [a,aint]=regress(y,X,0.02);a',aint'

58 >> yhat=y+sqrt(0.5)*randn(120,1);
>> [a,aint]=regress(yhat,X,0.02); >> a',aint‘ ％ a=[ ]' ans =

59 >> errorbar(1:6, a, aint(:,1)-a, aint(:,2)-a)
>> yhat=y+sqrt(0.1)*randn(120,1); >>[a,aint]=regress(yhat,X,0.02);

60 8.3.3 非线性函数的最小二乘参数 估计与区间估计 r为参数下的残 差构成的向量。 J为各个Jacobi 行向量构成的 矩阵。

61 例： >> f=inline('a(1)*exp(-a(2)*x)+a(3)*exp(-a(4)*x).*sin(a(5)*x)','a','x'); >> x=0:0.1:10; y=f([0.12,0.213,0.54,0.17,1.23],x); >> [a,r,j]=nlinfit(x,y,f,[1;1;1;1;1]); a

62 ans = >> ci=nlparci(a,r,j) ％ [0.12,0.213,0.54,0.17,1.23] ci =

63 >> y=f([0.12,0.213,0.54,0.17,1.23],x)+0.02*rand(size(x));
>> [a,r,j]=nlinfit(x,y,f,[1;1;1;1;1]); a' ans = >> ci=nlparci(a,r,j) ci = >> errorbar(1:5,a,ci(:,1)-a,ci(:,2)-a)

64 例： >> a=[1;1;1;1;1;1]';
>> f=inline(['(a(1)*x(:,1).^3+a(2)).*sin(a(3)*x(:,2) ' ，... '.*x(:,3))+(a(4)*x(:,3).^3+a(5)*x(:,3)+a(6))'],'a','x'); >> X=randn(120,3); y=f(a,X)+sqrt(0.2)*randn(120,1); >> [ahat,r,j]=nlinfit(X,y,f,[0;2;3;2;1;2]); ahat ahat =

65 >> ci=nlparci(ahat,r,j); ci ％置信区间
>> errorbar(1:6,ahat,ci(:,1)-ahat,ci(:,2)-ahat) >> y1=f(ahat,X);plot([y y1]) ％绘制曲线

66 8.4 统计假设检验 8.4.1 正态分布的均值假设检验 若未知正态分布的标准差时，可用此函数。
H为假设检验的结论，当H＝0时表示不拒绝H0假设，否则表示拒绝该假设。 s为接受假设的概率值， 为其均值的置信区间。 若未知正态分布的标准差时，可用此函数。

67 例：设某车间用一台包装机包装葡萄糖，包得的袋装糖重量是一个随机数，它服从正态分布。当机器正常时，其均值为0. 5公斤，标准差为0
例：设某车间用一台包装机包装葡萄糖，包得的袋装糖重量是一个随机数，它服从正态分布。当机器正常时，其均值为0.5公斤，标准差为0.015。某日开工后检验包装机是否正常，随机地抽取它所包装的的糖9袋，称得净重为（公斤）0.497, 0.506, 0.518, 0.524, 0.498, 0.511, 0.52, 0.515, 0.512,问机器是否正常？ 解： （分析）总体均值、标准差已知，则可设样本的标准差为0.015，于是 问题就化为根据样本值来判断 还是 。为此提出假设：

68 >> x=[0.497, 0.506, 0.518, 0.524, 0.498, 0.511, 0.52, 0.515, 0.512]; >> [H,p,ci]=ztest(x,0.5,0.015,0.05) H = 1 p = %样本观察值的概率 ci = %置信区间，均值0.5在此区间之外 结果H＝1，说明在0.05的水平下，拒绝原假设，即认为这天包装机工作不正常。

69 例：某种电子元件的寿命X（以小时计）服从正态分布，均值、方差均未知。现测得16只元件的寿命如下：
, 问是否有理由认为元件的平均寿命大于225（小时）： 解：按题意需做如下假设： 取

70 >> x=[ ]; >> [H,p,ci]=ttest(x,225,0.05) H = p = 0.6677 ci = %均值225在该置信区间内 结果表明，H＝0，即在显著水平为0.05的情况下，不能拒绝原假设。即认为元件的平均寿命不大于225小时。

71 8.4.2 正态分布假设检验 由随机样本判定分布是否为正态分布，可用下面两个假设算法的函数。
s为接受假设的概率值，s越接近于0，则可以拒绝是正态分布的原假设.

72 >> [H,p]=jbtest(X,0.05) %P为接受假设的概率值，P越接近于0，则可以拒绝是正态分布的原假设； H =
例： >> X=[216,203,197,208,206,209,206,208,202,203,206,213,218,207,208,... 202,194,203,213,211,193,213,208,208,204,206,204,206,208,209,... 213,203,206,207,196,201,208,207,213,208,210,208,211,211,214,... 220,211,203,216,224,211,209,218,214,219,211,208,221,211,218,... 218,190,219,211,208,199,214,207,207,214,206,217,214,201,212,... 213,211,212,216,206,210,216,204,221,208,209,214,214,199,204,... 211,201,216,211,209,208,209,202,211,207,202,205,206,216,206,... 213,206,207,200,198,200,202,203,208,216,206,222,213,209,219]; >> [H,p]=jbtest(X,0.05) %P为接受假设的概率值，P越接近于0，则可以拒绝是正态分布的原假设； H = p = 0.7281

73 >> [mu1,sig1,mu_ci,sig_ci]=normfit(X,0.05); mu=[mu1,mu_ci']
％该分布的均值及置信区间 >> sig=[sig1, sig_ci'] sig = ％该分布的方差及置信区间

74 %P为接受假设的概率值，P越接近于0，则可以拒绝是正态分布的原假设；c为测试统计量的值，d为是否拒绝原假设的临界值,c>d, 故拒绝。
例： >> r=gamrnd(1,3,400,1); [H,p,c,d]=jbtest(r,0.05) H = 1 p = c = d = 5.9915 %P为接受假设的概率值，P越接近于0，则可以拒绝是正态分布的原假设；c为测试统计量的值，d为是否拒绝原假设的临界值,c>d, 故拒绝。

75 8.4.3 其它分布的Kolmogorov-Smirnov 检验
其中cdffun为两列的值，第一列为自变量，第二列为对应的分布函数的值。

76 例： >> r=gamrnd(1,3,400,1); alam=gamfit(r) alam = 检验： >> r=sort(r); >> [H0,p]=kstest(r,[r gamcdf(r,alam(1),alam(2))],0.05) H0 = p = 0.6067

77 8.5方差分析及计算机求解 8.5.1 单因子方差分析 对一些观察来说，只有一个外界因素可能对观测的现象产生影响。
8.5方差分析及计算机求解 单因子方差分析 对一些观察来说，只有一个外界因素可能对观测的现象产生影响。 单因素方差分析是比较两组或多组数据的均值，它返回原假设—均值相等的概率,若p值接近于0，则原假设受到怀疑，说明至少有一列均值与其余列均值有明显不同。 X为需要分析的数据，每一列对应于随机分配的一个组的测试数据，这样会返回概率p，tab为方差分析表 。stats为统计结果量，为结构变量，包括每组均值等。

78 单因子方差分析表

79 例：

80 建立A矩阵，并求各列的均值。 >> A=[5,4,6,7,9; 8,6,4,4,3; 7,6,4,6,5; 7,3,5,6,7; 10,5,4,3,7; 8,6,3,5,6]; >> mean(A) ans = >> [p,tbl,stats]=anova1(A) %单因子方差分析 p = %<0.02或0.05，应拒绝给出的假设，有影响。 tbl = 'Source' 'SS' 'df' 'MS' 'F' 'Prob>F' 'Columns' [ ] [ 4] [9.1167] [3.8960] [0.0136] 'Error' [ ] [25] [2.3400] [] [] 'Total' [ ] [29] [] [] []

81 stats = gnames: [5x1 char] n: [ ] source: 'anova1' means: [ ] df: 25 s: 单因子方差表 盒式图

82 例：设有3台机器，用来生产规格相同的铝合金薄板。取样测量薄板的厚度，精确至‰厘米。得结果如下：
机器1： 机器2： 机器3： 检验各台机器所生产的薄板的厚度有无显著的差异？ >> X=[ ; ； ]; >> P=anova1(X') P = 1.3431e-005

83 8.5.2 双因子方差分析 如果有两种因子可能影响到某现象的统计规律，则应该引入双因子方差分析的概念。这时观察值可表示为一个三维数组。
根据双因子的特点，可以引入3个假设

84

85 双因素方差表

86 表中记号的定义 求解双因子方差分析问题:

87 例：比较 3 种松树在4 个不同地区的生长情况有无差别,在每个地区对每种松树随机地选择 5 株，测量它们的胸径，对它们进行双因子方差分析。

88 >> anova2(B‘,5); ％5表示每一单元观察点的数目
28,22,25,19,26,30,26,26,20,28,19,24,19,25,29,17,21,18,26,23; 18,10,12,22,13,15,21,22,14,12,23,25,19,13,22,16,12,23,22,19]; >> anova2(B‘,5); ％5表示每一单元观察点的数目 小， 很大，所以没有理由拒绝另外两个假设。故得出结论：地区对树的胸径无显著影响，不同区域对不同树种的胸径观测结果也无显著影响。

89 计算均值： >> C=[]; >> for i=1:3 for j=1:4 C(i,j)=mean(B(i,[1:5]+(j-1)*5)); end, end >> C=[C; mean(C)]; C=[C mean(C')'] C =

90 8.5.3 多因子方差分析

91

92 >> syms x y=dsolve('D2y-(2-1/x)*Dy+(1-1/x)*y=x^2*exp(-5*x)',... 'y(1)=pi','y(pi)=1','x') >> vpa(y,10) ans = e-3*exp(-5.*x)*(6.*ei(1,6.*x)*exp(6.*x) *x+36.*x^2) *exp(x) *log(x)*exp(x) >> x1=0.5:0.01:4; y1=subs(y,x,x1); >> plot(x1,y1,1,pi,'o',pi,1,'o')