数学实验概率论与数理统计问题的求解省公共课一等奖全国赛课获奖课件.pptx

第六讲 概率论与数理统计问题求解（下）概率分布与伪随机数生成统计量分析数理统计分析方法及计算机实现统计假设检验方差分析及计算机求解第1页 8.1概率分布与伪随机数生成 8.1.1 概率密度函数与分布函数概述第2页通用函数计算概率密度函数值 函数 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)第3页例:计算正态分布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第4页 随机变量累积概率值(分布函数值)通用函数cdf用来计算随机变量概率之和（累积概率值）函数 cdf格式 cdf(name，K，A)cdf(name，K，A，B)cdf(name，K，A，B，C)说明 返回以name为分布、随机变量XK概率之和累积概率值，name为分布函数名.第5页例:求标准正态分布随机变量X落在区间(-，0.4)内概率。解：cdf(norm,0.4,0,1)ans=0.6554例：求自由度为16卡方分布随机变量落在0，6.91内概率。解：cdf(chi2,6.91,16)ans=0.0250第6页随机变量逆累积分布函数 MATLAB中逆累积分布函数是已知，求x。命令 icdf 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)第7页例:在标准正态分布表中，若已知F=0.6554,求X解：icdf(norm,0.6554,0,1)ans=0.3999例：公共汽车门高度是按成年男子与车门顶碰头机会不超出1%设计。设男子身高X（单位：cm）服从正态分布N（175，6），求车门最低高度。解：设h为车门高度，X为身高。求满足条件 FXh=0.99，即 FX=0.01故 h=icdf(norm,0.99,175,6)h=188.9581第8页8.1.2 常见分布概率密度函数与分布函数 8.1.2.1 Poisson分布其要求x是正整数。第9页其中：x为选定一组横坐标向量，y为x各点处概率密度函数值。第10页例：绘制 l 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)第11页8.1.2.2 正态分布正态分布概率密度函数为:第12页例：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)第13页8.1.2.3 分布第14页例：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)第15页8.1.2.4 分布（卡方分布）其为一特殊 分布，a=k/2,l l=1/2。第16页例：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)第17页8.1.2.5 分布概率密度函数为:其为参数k函数，且k为正整数。第18页例：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)第19页8.1.2.6 Rayleigh分布第20页例：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)第21页8.1.2.7 F 分布其为参数p,q函数，且p,q均为正整数。第22页例：分别绘制（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=1 2 3 3 4;q1=1 1 1 2 1;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)第23页8.1.3 概率问题求解图4-9第24页例：b=1;p1=raylcdf(0.2,b);p2=raylcdf(2,b);P1=p2-p1P1=0.8449 p1=raylcdf(1,b);P2=1-p1P2=0.6065第25页例：syms x y;f=x2+x*y/3;P=int(int(f,x,0,1/2),y,0,1/2)P=5/192 syms x y;f=x2+x*y/3;P=int(int(f,x,0,1),y,0,2)P=1第26页8.1.4 随机数与伪随机数第27页第28页例：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)第29页8.2 统计量分析 8.2.1 随机变量均值与方差第30页例：均值 syms x;syms a lam positive p=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第31页已知一组随机变量样本数据组成向量:求该向量各个元素均值、方差和标准差、中位数medianmedian第32页例：生成一组 30000 个正态分布随机数，使其均值为 0.5，标准差为1.5，分析数据实际均值、方差和标准差，假如减小随机变量个数，会有什么结果？p=normrnd(0.5,1.5,30000,1);mean(p),var(p),std(p)ans=0.4879 2.2748 1.5083300个随机数 p=normrnd(0.5,1.5,300,1);mean(p),var(p),std(p)ans=0.4745 1.9118 1.3827可见在进行较准确统计分析时不能选择太小样本点。第33页例：m,s=raylstat(0.45)m=0.5640s=0.0869第34页8.2.2 随机变量矩第35页例：求解原点矩 syms x;syms a lam positive;p=lama*x(a-1)/gamma(a)*exp(-lam*x);for n=1:5,m=int(xn*p,x,0,inf),endm=1/lam*a m=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)有规律第36页 syms n;m=simple(int(x)n*p,x,0,inf)直接求出m=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=0s=a/lam2 s=2*a/lam3s=3*a*(a+2)/lam4s=4*a*(5*a+6)/lam5s=5*a*(3*a2+26*a+24)/lam6 好像无规律第37页第38页例：考虑前面随机数，能够用下面语句得出随机数各阶矩。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,BA=0.5066 2.4972 3.5562 18.7530 41.5506B=0 2.2405 0.0212 15.1944 0.0643第39页求各阶距理论值：syms x;A1=;B1=;p=1/(sqrt(2*pi)*1.5)*exp(-(x-0.5)2/(2*1.52);for i=1:5 A1=A1,vpa(int(xi*p,x,-inf,inf),12);B1=B1,vpa(int(x-0.5)i*p,x,-inf,inf),12);end A1,B1A1=.500000000001,2.50000000000,3.50000000001,18.6250000000,40.8125000000 B1=0,2.25000000000,0,15.1875000000,0第40页8.2.3 多变量随机数协方差分析第41页第42页例：p=randn(30000,4);cov(p)ans=1.0033 0.0131 0.0036 0.0020 0.0131 1.0110 0.0061 -0.0154 0.0036 0.0061 1.0055 -0.0004 0.0020 -0.0154 -0.0004 0.9881第43页8.2.4 多变量正态分布联合概率密度即分布函数第44页例：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*161*41)surf(X,Y,P)第45页 对协方差矩阵进行处理，可计算出新联合概率密度函数。Sigma2=diag(diag(Sigma2);%消除协方差矩阵非对角元素 p=mvnpdf(xy,mu1,Sigma2);P=reshape(p,size(X);surf(X,Y,P)R为m行n列。第46页例：mu1=-1,2;Sigma2=1 1;1 3;R1=mvnrnd(mu1,Sigma2,);plot(R1(:,1),R1(:,2),o)Sigma2=diag(diag(Sigma2);figure;R2=mvnrnd(mu1,Sigma2,);plot(R2(:,1),R2(:,2),o)第47页8.3数理统计分析方法及计算机实现 8.3.1 参数预计与区间预计 不论总体X分布函数F（x；）类型已知或未知，我们总是需要去预计一些未知参数或数字特征，这就是参数预计问题.即参数预计就是从样本（X1，X2，Xn）出发，结构一些统计量 X1，X2，Xn）（i=1，2，k）去预计总体X中一些参数（或数字特征）（i=1，2，k）.这么统计量称为预计量预计量.第48页1、点预计、点预计：结构（X1，X2，Xn）函数 (X1，X2，Xn）作为参数 点预计量，称统计量 为总体X参数 点预计量.2.区间预计区间预计：结构两个函数 (X1，X2，Xn）和 (X1，X2，Xn）做成区间，把这 （）作为参数 区间预计.第49页区间预计求法区间预计求法 设总体X分布中含有未知参数 ，若对于给定概率 ，存在两个统计量 (X1，X2，Xn）和 (X1，X2，Xn）,使得 则称随机区间 为参数 置信水平为 置信区间置信区间,称 为置信下限置信下限,称 为置信上限置信上限.第50页 由极大拟然法预计出该分布均值、方差 及其置信区间。置信度越大，得出置信区间越小，即得出结果越靠近于真值。还有gamfit(),raylfit(),poissfit()，unifit()（均匀分布）等参数预计函数第51页例：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 AA=0.9000 1.5137 1.5123 1.5152 2.9825 2.9791 2.9858 0.9200 1.5137 1.5126 1.5149 2.9825 2.9798 2.9851 0.9500 1.5137 1.5130 1.5144 2.9825 2.9808 2.9841 0.9800 1.5137 1.5135 1.5140 2.9825 2.9818 2.9831第52页 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=1.4795 1.4725 1.4865 2.9129 2.8960 2.9299 1.4218 1.4198 1.4238 3.1676 3.1623 3.1729 1.4898 1.4891 1.4904 3.0425 3.0409 3.0442 1.4998 1.4996 1.5000 3.0054 3.0049 3.0059 1.5006 1.5005 1.5007 2.9968 2.9966 2.9969 要到达参数预计效果良好，随机数不能选得太少，也不能选得太多，此例中为30000为好。第53页8.3.2 多元线性回归与区间预计第54页第55页第56页例：a=1-1.232 2.23 2 4 3.792;X=randn(120,6);y=X*a;a1=inv(X*X)*X*y;a1ans=1.0000 -1.2320 2.2300 2.0000 4.0000 3.7920 a,aint=regress(y,X,0.02);a,aintans=1.0000 -1.2320 2.2300 2.0000 4.0000 3.7920ans=1.0000 -1.2320 2.2300 2.0000 4.0000 3.7920 1.0000 -1.2320 2.2300 2.0000 4.0000 3.7920第57页 yhat=y+sqrt(0.5)*randn(120,1);a,aint=regress(yhat,X,0.02);a,aint a=1-1.232 2.23 2 4 3.792ans=1.0576 -1.3280 2.1832 2.0151 4.0531 3.7749ans=0.8800 -1.5107 2.0284 1.8544 3.8788 3.6221 1.2353 -1.1453 2.3379 2.1757 4.2274 3.9276第58页 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)第59页8.3.3 非线性函数最小二乘参数预计与区间预计r为参数下残差组成向量。J为各个Jacobi行向量组成矩阵。第60页例：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第61页ans=0.11999999763418 0.21299999458274 0.54000000196818 0.17000000068705 1.22999999996315 ci=nlparci(a,r,j)0.12,0.213,0.54,0.17,1.23ci=0.11999999712512 0.11999999814323 0.21299999340801 0.21299999575747 0.54000000124534 0.54000000269101 0.17000000036077 0.17000000101332 1.22999999978603 1.23000000014028第62页 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.12655784086874 0.17576593556541 0.54363873794463 0.17129712329146 1.23139632101927 ci=nlparci(a,r,j)ci=0.12240417108574 0.130711510651740.16754837168468 0.183983499446140.53737093469422 0.549906541195040.16845014477426 0.174144101808661.22983289563708 1.23295974640145 errorbar(1:5,a,ci(:,1)-a,ci(:,2)-a)第63页例：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.99166464884539 1.04776526972943 0.97668595800756 1.02022345889541 0.88639528713563 1.09317291667891第64页 ci=nlparci(ahat,r,j);ci 置信区间ci=0.89133624667624 1.09199305101455 0.86664749663205 1.22888304282680 0.83628948119418 1.11708243482094 0.98466523279168 1.05578168499914 0.73055684224143 1.04223373202984 0.99932407018303 1.18702176317478 errorbar(1:6,ahat,ci(:,1)-ahat,ci(:,2)-ahat)y1=f(ahat,X);plot(y y1)绘制曲线第65页8.4 统计假设检验8.4.1 正态分布均值假设检验 H为假设检验结论，当H0时表示不拒绝H0假设，不然表示拒绝该假设。s为接收假设概率值，为其均值置信区间。若未知正态分布标准差时，可用此函数。第66页例：设某车间用一台包装机包装葡萄糖，包得袋装糖重量是一个随机数，它服从正态分布。当机器正常时，其均值为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，于是 问题就化为依据样本值来判断 还是 。为此提出假设：第67页 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.5014 0.5210%置信区间，均值0.5在此区间之外 结果H1，说明在0.05水平下，拒绝原假设，即认为这天包装机工作不正常。第68页例：某种电子元件寿命X（以小时计）服从正态分布，均值、方差均未知。现测得16只元件寿命以下：159 280 101 212 224 379 179 264 222 262 168 250 149 260 485 170,问是否有理由认为元件平均寿命大于225（小时）：解：按题意需做以下假设：取第69页 x=159 280 101 212 224 379 179 264 222 262 168 250 149 260 485 170;H,p,ci=ttest(x,225,0.05)H=0p=0.6677ci=185.3622 285.1378%均值225在该置信区间内 结果表明，H0，即在显著水平为0.05情况下，不能拒绝原假设。即认为元件平均寿命小于225小时。第70页8.4.2 正态分布假设检验 由随机样本判定分布是否为正态分布，可用下面两个假设算法函数。s为接收假设概率值，s越靠近于0，则能够拒绝是正态分布原假设.第71页例：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=0 p=0.7281第72页 mu1,sig1,mu_ci,sig_ci=normfit(X,0.05);mu=mu1,mu_cimu=208.8167 207.6737 209.9596该分布均值及置信区间 sig=sig1,sig_cisig=6.3232 5.6118 7.2428该分布方差及置信区间第73页例：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,故拒绝。第74页8.4.3 其它分布Kolmogorov-Smirnov 检验 此函数（Kolmogorov-Smirnov 算法）可对任意已知分布函数进行有效假设检验。其中cdffun为两列值，第一列为自变量，第二列为对应分布函数值。第75页例：r=gamrnd(1,3,400,1);alam=gamfit(r)alam=0.9708 3.1513检验：r=sort(r);H0,p=kstest(r,r gamcdf(r,alam(1),alam(2),0.05)H0=0p=0.6067第76页 8.5方差分析及计算机求解 8.5.1 单因子方差分析 对一些观察来说，只有一个外界原因可能对观察现象产生影响。单原因方差分析是比较两组或多组数据均值，它返回原假设均值相等概率,若p值靠近于0，则原假设受到怀疑，说明最少有一列均值与其余列均值有显著不一样。X为需要分析数据，每一列对应于随机分配一个组测试数据，这么会返回概率p，tab为方差分析表。stats为统计结果量，为结构变量，包含每组均值等。第77页单因子方差分析表第78页例：第79页建立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.5000 5.0000 4.3333 5.1667 6.1667 p,tbl,stats=anova1(A)%单因子方差分析p=0.0136%F Columns 36.4667 4 9.1167 3.8960 0.0136 Error 58.5000 25 2.3400 Total 94.9667 29 第80页stats=gnames:5x1 char n:6 6 6 6 6 source:anova1 means:7.5000 5 4.3333 5.1667 6.1667 df:25 s:1.5297单因子方差表 盒式图第81页例：设有3台机器，用来生产规格相同铝合金薄板。取样测量薄板厚度，准确至厘米。得结果以下：机器1：0.236 0.238 0.248 0.245 0.243机器2：0.257 0.253 0.255 0.254 0.261机器3：0.258 0.264 0.259 0.267 0.262检验各台机器所生产薄板厚度有没有显著差异？X=0.236 0.238 0.248 0.245 0.243;0.257 0.253 0.255 0.254 0.261；0.258 0.264 0.259 0.267 0.262;P=anova1(X)P=1.3431e-005第82页8.5.2 双因子方差分析 假如有两种因子可能影响到某现象统计规律，则应该引入双因子方差分析概念。这时观察值可表示为一个三维数组。依据双因子特点，能够引入3个假设第83页第84页双原因方差表第85页表中记号定义求解双因子方差分析问题:第86页例：比较 3 种松树在4 个不一样地域生长情况有没有差异,在每个地域对每种松树随机地选择 5 株，测量它们胸径，对它们进行双因子方差分析。第87页 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表示每一单元观察点数目 小(有影响)，很大(无影响)，所以没有理由拒绝另外两个假设。故得出结论：树之间有显著差异,地域对树胸径无显著影响，不一样区域对不一样树种胸径观察结果也无显著影响。第88页计算均值：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=19.6000 20.0000 21.0000 18.8000 19.8500 24.0000 26.0000 23.21.0000 23.5500 15.0000 16.8000 20.4000 18.4000 17.6500 19.5333 20.9333 21.5333 19.4000 20.3500第89页8.5.3 多因子方差分析第90页第91页 syms xy=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)第92页