数学实验6-2概率论与数理统计问题的求解.ppt

第六讲概率论与数理统计问题的求解（下）,概率分布与伪随机数生成统计量分析数理统计分析方法及计算机实现统计假设检验方差分析及计算机求解,8.1概率分布与伪随机数生成8.1.1概率密度函数与分布函数概述,通用函数计算概率密度函数值,函数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=PX=K=pdf(bino，K，n，p),例:计算正态分布N（0，1）的随机变量X在点0.6578的密度函数值。解：pdf(norm,0.6578,0,1)ans=0.3213例:自由度为8的卡方分布，在点2.18处的密度函数值。解：pdf(chi2,2.18,8)ans=0.0363,随机变量的累积概率值(分布函数值),通用函数cdf用来计算随机变量的概率之和（累积概率值）函数cdf格式cdf(name，K，A)cdf(name，K，A，B)cdf(name，K，A，B，C)说明返回以name为分布、随机变量XK的概率之和的累积概率值，name为分布函数名.,例:求标准正态分布随机变量X落在区间(-，0.4)内的概率。解：cdf(norm,0.4,0,1)ans=0.6554例：求自由度为16的卡方分布随机变量落在0，6.91内的概率。解：cdf(chi2,6.91,16)ans=0.0250,随机变量的逆累积分布函数,MATLAB中的逆累积分布函数是已知，求x。命令icdf计算逆累积分布函数格式icdf(name，K，A)icdf(name，K，A，B)icdf(name，K，A，B，C)说明返回分布为name，参数为a1,a2,a3，累积概率值为P的临界值，这里name与前面相同。如果F=cdf(name，X，A，B，C)，则X=icdf(name，F，A，B，C),例:在标准正态分布表中，若已知F=0.6554,求X解：icdf(norm,0.6554,0,1)ans=0.3999例：公共汽车门的高度是按成年男子与车门顶碰头的机会不超过1%设计的。设男子身高X（单位：cm）服从正态分布N（175，6），求车门的最低高度。解：设h为车门高度，X为身高。求满足条件FXh=0.01故h=icdf(norm,0.99,175,6)h=188.9581,8.1.2常见分布的概率密度函数与分布函数8.1.2.1Poisson分布,其要求x是正整数。,其中：x为选定的一组横坐标向量，y为x各点处的概率密度函数值。,例：绘制l=1,2,5,10时Poisson分布的概率密度函数与概率分布函数曲线。x=0:15;y1=;y2=;lam1=1,2,5,10;fori=1:length(lam1)y1=y1,poisspdf(x,lam1(i);y2=y2,poisscdf(x,lam1(i);endplot(x,y1),figure;plot(x,y2),8.1.2.2正态分布,正态分布的概率密度函数为:,例：x=-5:.02:5;y1=;y2=;mu1=-1,0,0,0,1;sig1=1,0.1,1,10,1;sig1=sqrt(sig1);fori=1:length(mu1)y1=y1,normpdf(x,mu1(i),sig1(i);y2=y2,normcdf(x,mu1(i),sig1(i);endplot(x,y1),figure;plot(x,y2),8.1.2.3分布,例：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;fori=1:length(a1)y1=y1,gampdf(x,a1(i),lam1(i);y2=y2,gamcdf(x,a1(i),lam1(i);endplot(x,y1),figure;plot(x,y2),8.1.2.4分布（卡方分布）,其为一特殊的分布，a=k/2,l=1/2。,例：x=-eps:-0.02:-0.5,0:0.02:2;x=sort(x);k1=1,2,3,4,5;y1=;y2=;fori=1:length(k1)y1=y1,chi2pdf(x,k1(i);y2=y2,chi2cdf(x,k1(i);endplot(x,y1),figure;plot(x,y2),8.1.2.5分布,概率密度函数为:,其为参数k的函数，且k为正整数。,例：x=-5:0.02:5;k1=1,2,5,10;y1=;y2=;fori=1:length(k1)y1=y1,tpdf(x,k1(i);y2=y2,tcdf(x,k1(i);endplot(x,y1),figure;plot(x,y2),8.1.2.6Rayleigh分布,例：x=-eps:-0.02:-0.5,0:0.02:5;x=sort(x);b1=.5,1,3,5;y1=;y2=;fori=1:length(b1)y1=y1,raylpdf(x,b1(i);y2=y2,raylcdf(x,b1(i);endplot(x,y1),figure;plot(x,y2),8.1.2.7F分布,其为参数p,q的函数，且p,q均为正整数。,例：分别绘制（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=12334;q1=11121;y1=;y2=;fori=1:length(p1)y1=y1,fpdf(x,p1(i),q1(i);y2=y2,fcdf(x,p1(i),q1(i);endplot(x,y1),figure;plot(x,y2),8.1.3概率问题的求解,图4-9,例：b=1;p1=raylcdf(0.2,b);p2=raylcdf(2,b);P1=p2-p1P1=0.8449p1=raylcdf(1,b);P2=1-p1P2=0.6065,例：symsxy;f=x2+x*y/3;P=int(int(f,x,0,1/2),y,0,1/2)P=5/192symsxy;f=x2+x*y/3;P=int(int(f,x,0,1),y,0,2)P=1,8.1.4随机数与伪随机数,例：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),8.2统计量分析8.2.1随机变量的均值与方差,例：均值symsx;symsalampositivep=lama*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/lam2,已知一组随机变量样本数据构成的向量:,求该向量各个元素的均值、方差和标准差、中位数median,例：生成一组30000个正态分布随机数，使其均值为0.5，标准差为1.5，分析数据实际的均值、方差和标准差，如果减小随机变量个数，会有什么结果？p=normrnd(0.5,1.5,30000,1);mean(p),var(p),std(p)ans=0.48792.27481.5083300个随机数p=normrnd(0.5,1.5,300,1);mean(p),var(p),std(p)ans=0.47451.91181.3827可见在进行较精确的统计分析时不能选择太小的样本点。,例：m,s=raylstat(0.45)m=0.5640s=0.0869,8.2.2随机变量的矩,例：求解原点矩symsx;symsalampositive;p=lama*x(a-1)/gamma(a)*exp(-lam*x);forn=1:5,m=int(xn*p,x,0,inf),endm=1/lam*am=1/lam2*a*(a+1)m=1/lam3*a*(a+1)*(a+2)m=1/lam4*a*(a+1)*(a+2)*(a+3)m=1/lam5*a*(a+1)*(a+2)*(a+3)*(a+4)有规律,symsn;m=simple(int(x)n*p,x,0,inf)直接求出m=lam(-n)*gamma(n+a)/gamma(a)forn=1:6,s=simple(int(x-1/lam*a)n*p,x,0,inf),end中心距s=0s=a/lam2s=2*a/lam3s=3*a*(a+2)/lam4s=4*a*(5*a+6)/lam5s=5*a*(3*a2+26*a+24)/lam6好像无规律,例：考虑前面的随机数，可以用下面的语句得出随机数的各阶矩。A=;B=;p=normrnd(0.5,1.5,30000,1);n=1:5;forr=n,A=A,sum(p.r)/length(p);B=B,moment(p,r);endA,BA=0.50662.49723.556218.753041.5506B=02.24050.021215.19440.0643,求各阶距的理论值：symsx;A1=;B1=;p=1/(sqrt(2*pi)*1.5)*exp(-(x-0.5)2/(2*1.52);fori=1:5A1=A1,vpa(int(xi*p,x,-inf,inf),12);B1=B1,vpa(int(x-0.5)i*p,x,-inf,inf),12);endA1,B1A1=.500000000001,2.50000000000,3.50000000001,18.6250000000,40.8125000000B1=0,2.25000000000,0,15.1875000000,0,8.2.3多变量随机数的协方差分析,例：p=randn(30000,4);cov(p)ans=1.00330.01310.00360.00200.01311.01100.0061-0.01540.00360.00611.0055-0.00040.0020-0.0154-0.00040.9881,8.2.4多变量正态分布的联合概率密度即分布函数,例：mu1=-1,2;Sigma2=11;13;%输入均值向量和协方差矩阵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);Changesize(2501*161*41)surf(X,Y,P),对协方差矩阵进行处理，可计算出新的联合概率密度函数。Sigma2=diag(diag(Sigma2);%消除协方差矩阵的非对角元素p=mvnpdf(xy,mu1,Sigma2);P=reshape(p,size(X);surf(X,Y,P),R为m行n列。,例：mu1=-1,2;Sigma2=11;13;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),8.3数理统计分析方法及计算机实现8.3.1参数估计与区间估计,无论总体X的分布函数F（x；）的类型已知或未知，我们总是需要去估计某些未知参数或数字特征，这就是参数估计问题.即参数估计就是从样本（X1，X2，Xn）出发，构造一些统计量X1，X2，Xn）（i=1，2，k）去估计总体X中的某些参数（或数字特征）（i=1，2，k）.这样的统计量称为估计量.,1、点估计：构造（X1，X2，Xn）的函数(X1，X2，Xn）作为参数的点估计量，称统计量为总体X参数的点估计量.2.区间估计：构造两个函数(X1，X2，Xn）和(X1，X2，Xn）做成区间，把这（）作为参数的区间估计.,区间估计的求法,设总体X的分布中含有未知参数，若对于给定的概率，存在两个统计量(X1，X2，Xn）和(X1，X2，Xn）,使得则称随机区间为参数的置信水平为的置信区间,称为置信下限,称为置信上限.,由极大拟然法估计出该分布的均值、方差及其置信区间。置信度越大，得出的置信区间越小，即得出的结果越接近于真值。还有gamfit(),raylfit(),poissfit()，unifit()（均匀分布）等参数估计函数,例：p=gamrnd(1.5,3,30000,1);Pv=0.9,0.92,0.95,0.98;A=;fori=1:length(Pv)a,b=gamfit(p,Pv(i);A=A;Pv(i),a(1),b(:,1),a(2),b(:,2)endAA=0.90001.51371.51231.51522.98252.97912.98580.92001.51371.51261.51492.98252.97982.98510.95001.51371.51301.51442.98252.98082.98410.98001.51371.51351.51402.98252.98182.9831,num=300,3000,30000,300000,3000000;A=;fori=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);endA(:,2,3,4,5,6,7)ans=1.47951.47251.48652.91292.89602.92991.42181.41981.42383.16763.16233.17291.48981.48911.49043.04253.04093.04421.49981.49961.50003.00543.00493.00591.50061.50051.50072.99682.99662.9969要达到参数估计效果良好，随机数不能选得太少，也不能选得太多，此例中为30000为好。,8.3.2多元线性回归与区间估计,例：a=1-1.2322.23243.792;X=randn(120,6);y=X*a;a1=inv(X*X)*X*y;a1ans=1.0000-1.23202.23002.00004.00003.7920a,aint=regress(y,X,0.02);a,aintans=1.0000-1.23202.23002.00004.00003.7920ans=1.0000-1.23202.23002.00004.00003.79201.0000-1.23202.23002.00004.00003.7920,yhat=y+sqrt(0.5)*randn(120,1);a,aint=regress(yhat,X,0.02);a,ainta=1-1.2322.23243.792ans=1.0576-1.32802.18322.01514.05313.7749ans=0.8800-1.51072.02841.85443.87883.62211.2353-1.14532.33792.17574.22743.9276,errorbar(1:6,a,aint(:,1)-a,aint(:,2)-a)errorbar()用图形绘制参数估计的置信区间。yhat=y+sqrt(0.1)*randn(120,1);a,aint=regress(yhat,X,0.02);errorbar(1:6,a,aint(:,1)-a,aint(:,2)-a),8.3.3非线性函数的最小二乘参数估计与区间估计,r为参数下的残差构成的向量。J为各个Jacobi行向量构成的矩阵。,例：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,ans=0.119999997634180.212999994582740.540000001968180.170000000687051.22999999996315ci=nlparci(a,r,j)0.12,0.213,0.54,0.17,1.23ci=0.119999997125120.119999998143230.212999993408010.212999995757470.540000001245340.540000002691010.170000000360770.170000001013321.229999999786031.23000000014028,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);aans=0.126557840868740.175765935565410.543638737944630.171297123291461.23139632101927ci=nlparci(a,r,j)ci=0.122404171085740.130711510651740.167548371684680.183983499446140.537370934694220.549906541195040.168450144774260.174144101808661.229832895637081.23295974640145errorbar(1:5,a,ci(:,1)-a,ci(:,2)-a),例：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);ahatahat=0.991664648845391.047765269729430.976685958007561.020223458895410.886395287135631.09317291667891,ci=nlparci(ahat,r,j);ci置信区间ci=0.891336246676241.091993051014550.866647496632051.228883042826800.836289481194181.117082434820940.984665232791681.055781684999140.730556842241431.042233732029840.999324070183031.18702176317478errorbar(1:6,ahat,ci(:,1)-ahat,ci(:,2)-ahat)y1=f(ahat,X);plot(yy1)绘制曲线,8.4统计假设检验8.4.1正态分布的均值假设检验,H为假设检验的结论，当H0时表示不拒绝H0假设，否则表示拒绝该假设。s为接受假设的概率值，为其均值的置信区间。若未知正态分布的标准差时，可用此函数。,例：设某车间用一台包装机包装葡萄糖，包得的袋装糖重量是一个随机数，它服从正态分布。当机器正常时，其均值为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，于是问题就化为根据样本值来判断还是。为此提出假设：,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=1p=0.0248%样本观察值的概率ci=0.50140.5210%置信区间，均值0.5在此区间之外结果H1，说明在0.05的水平下，拒绝原假设，即认为这天包装机工作不正常。,例：某种电子元件的寿命X（以小时计）服从正态分布，均值、方差均未知。现测得16只元件的寿命如下：159280101212224379179264222262168250149260485170,问是否有理由认为元件的平均寿命大于225（小时）：解：按题意需做如下假设：取,x=159280101212224379179264222262168250149260485170;H,p,ci=ttest(x,225,0.05)H=0p=0.6677ci=185.3622285.1378%均值225在该置信区间内结果表明，H0，即在显著水平为0.05的情况下，不能拒绝原假设。即认为元件的平均寿命不大于225小时。,8.4.2正态分布假设检验,由随机样本判定分布是否为正态分布，可用下面两个假设算法的函数。s为接受假设的概率值，s越接近于0，则可以拒绝是正态分布的原假设.,例：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=0p=0.7281,mu1,sig1,mu_ci,sig_ci=normfit(X,0.05);mu=mu1,mu_cimu=208.8167207.6737209.9596该分布的均值及置信区间sig=sig1,sig_cisig=6.32325.61187.2428该分布的方差及置信区间,例：r=gamrnd(1,3,400,1);H,p,c,d=jbtest(r,0.05)H=1p=0c=504.2641d=5.9915%P为接受假设的概率值，P越接近于0，则可以拒绝是正态分布的原假设；c为测试统计量的值，d为是否拒绝原假设的临界值,cd,故拒绝。,8.4.3其它分布的Kolmogorov-Smirnov检验,此函数（Kolmogorov-Smirnov算法）可对任意已知分布函数进行有效的假设检验。其中cdffun为两列的值，第一列为自变量，第二列为对应的分布函数的值。,例：r=gamrnd(1,3,400,1);alam=gamfit(r)alam=0.97083.1513检验：r=sort(r);H0,p=kstest(r,rgamcdf(r,alam(1),alam(2),0.05)H0=0p=0.6067,8.5方差分析及计算机求解8.5.1单因子方差分析,对一些观察来说，只有一个外界因素可能对观测的现象产生影响。单因素方差分析是比较两组或多组数据的均值，它返回原假设均值相等的概率,若p值接近于0，则原假设受到怀疑，说明至少有一列均值与其余列均值有明显不同。X为需要分析的数据，每一列对应于随机分配的一个组的测试数据，这样会返回概率p，tab为方差分析表。stats为统计结果量，为结构变量，包括每组均值等。,单因子方差分析表,例：,建立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=7.50005.00004.33335.16676.1667p,tbl,stats=anova1(A)%单因子方差分析p=0.0136%FColumns36.466749.11673.89600.0136Error58.5000252.3400Total94.966729,stats=gnames:5x1charn:66666source:anova1means:7.500054.33335.16676.1667df:25s:1.5297单因子方差表盒式图,例：设有3台机器，用来生产规格相同的铝合金薄板。取样测量薄板的厚度，精确至厘米。得结果如下：机器1：0.2360.2380.2480.2450.243机器2：0.2570.2530.2550.2540.261机器3：0.2580.2640.2590.2670.262检验各台机器所生产的薄板的厚度有无显著的差异？X=0.2360.2380.2480.2450.243;0.2570.2530.2550.2540.261；0.2580.2640.2590.2670.262;P=anova1(X)P=1.3431e-005,8.5.2双因子方差分析,如果有两种因子可能影响到某现象的统计规律，则应该引入双因子方差分析的概念。这时观察值可表示为一个三维数组。根据双因子的特点，可以引入3个假设,双因素方差表,表中记号的定义,求解双因子方差分析问题:,例：比较3种松树在4个不同地区的生长情况有无差别,在每个地区对每种松树随机地选择5株，测量它们的胸径，对它们进行双因子方差分析。,B=23,15,26,13,21,25,20,21,16,18,21,17,16,24,27,14,17,19,20,24;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表示每一单元观察点的数目小(有影响)，很大(无影响)，所以没有理由拒绝另外两个假设。故得出结论：树之间有显著差异,地区对树的胸径无显著影响，不同区域对不同树种的胸径观测结果也无显著影响。,计算均值：C=;fori=1:3forj=1:4C(i,j)=mean(B(i,1:5+(j-1)*5);end,endC=C;mean(C);C=Cmean(C)C=19.600020.000021.000018.800019.850024.000026.000023.200021.000023.550015.000016.800020.400018.400017.650019.533320.933321.533319.400020.3500,8.5.3多因子方差分析,symsxy=dsolve(D2y-(2-1/x)*Dy+(1-1/x)*y=x2*exp(-5*x),.y(1)=pi,y(pi)=1,x)vpa(y,10)ans=.7716049383e-3*exp(-5.*x)*(6.*ei(1,6.*x)*exp(6.*x)+11.+30.*x+36.*x2)+1.155578411*exp(x)-.9717266135*log(x)*exp(x)x1=0.5:0.01:4;y1=subs(y,x,x1);plot(x1,y1,1,pi,o,pi,1,o),