matlab概率统计实验
matlab概率统计实验
9.1 实验（i）：galton钉板试验
9.1.1 实验与观察: galton钉板模型和二项分布
1. 动画模拟calton钉板试验
【 rand('seed',1), u=rand(1,6) 】
【 rand('seed',2),u=rand(1,6) 】
观察程序zxy9_1.m
【 clear，clf, m=100;n=5;y0=2; %设置参数。
ballnum=zeros(1,n+1);
p=0.5;q=1-p;
for i=n+1:-1:1 %创建钉子的坐标x,y 。
x(i,1)=0.5*(n-i+1);y(i,1)=(n-i+1)+y0;
for j=2:i
x(i,j)=x(i,1)+(j-1)*1;y(i,j)=y(i,1);
end
end
mm=moviein(m); % 动画开始，模拟小球下落路径。
for i=1:m
s=rand(1,n); %产生n个随机数。
xi=x(1,1);yi=y(1,1);k=1;l=1; % 小球遇到第一个钉子。
for j=1:n
plot(x(1:n,:),y(1:n,:),'o',x(n+1,:),y(n+1,:),'.-'), % 画钉子的位置 。
axis([-2 n+2 0 y0+n+1]),hold on
k=k+1; % 小球下落一格 。
if s(j)>p
l=l+0; %小球左移 。
else
l=l+1; %小球右移 。
end
xt=x(k,l);yt=y(k,l); %小球下落点的坐标。
h=plot([xi,xt],[yi,yt]);axis([-2 n+2 0 y0+n+1]) %画小球运动轨迹。
xi=xt;yi=yt;
end
ballnum(l)=ballnum(l)+1; %计数。
ballnum1=3*ballnum./m;
bar([0:n],ballnum1),axis([-2 n+2 0 y0+n+1]) %画各格子的频率 。
mm(i)=getframe; %存储动画数据。
hold off
end
movie(mm,1) %播放动画1次。 】
2. 用二项分布描述galton钉板模型
【 x = binornd(5,0.5,1,10) 】
3. 数学期望和平均收益
【 n=5;p=0.5; m=5000;
rand('seed',5); r = binornd(n,p,1,m); %模拟投球。
f=[5, 1, 0.2, 0.2, 1, 5]; %奖品的价值向量。
s=[];
for k=1:n+1 %计算第0～n格奖品价值。
u=[];
u=find(r==(k-1)); %计算落入第k-1格的小球下标，并存于向量u中。
s=[f(k)*length(u),s]; %计算相应的奖品价值，并存于向量s中。
end
mean_return=sum(s)/m %计算一次抽奖的平均回报 。 】
{ mean_return = 0.7506 }
【 f=[5, 1, 0.2, 0.2, 1, 5]; x=[0:1:n]; pu=binopdf(x,n,p);ef=sum(f.*pu) 】
9.1.2 应用、思考和练习
1. 二项分布的应用模型
◆电力供应问题 。
提示：
【 x=linspace(90,160,10); r = binornd(200,0.6,1,960); hist(r,x) ;
[n,x]=hist(r,x); 】
◆ 废品问题。
下面的程序可作为参考 。
【 clear,clf ,n=50;r=linspace(0,1,n);s=linspace(0,10,n);
[r,s]=meshgrid(r,s); x=[0 1 2 3 4 5];p = binopdf(x,5,0.5);
for k=1:n
f=[s(k,l) s(k,l)*r(k,l) s(k,l)*r(k,l)^2 s(k,l)*r(k,l)^2 s(k,l)*r(k,l) s(k,l)];
for l=1:n
z(k,l)=sum(f.*p);
if z(k,l)>=1 z(k,l)=nan; end
end
end
contour(r,s,z) 】
3. 单服务台定长服务时间排队系统的计算机模拟
4. 随机变量的模拟：逆概率法
图9.4 逆概率法模拟随机数原理。
【 clf
mu=0;sigma=1;a=mu-4*sigma;b=mu+4*sigma;
t=linspace(a,b,30);
f1=normcdf(t,mu,sigma);
for i=1:3
hold on
u=rand(1);
f=norminv(u,mu,sigma);
plot(t,f1,[0 0],[0,b]),hold on,
plot(0,u,'ro'),plot([0,f],[u,u]);plot([f,f],[u,0]);
plot(f,0,'ro'), axis([-4 4 0 1])
end 】
9.2 实验(ⅱ) ：热轧机的调整
9.2.1实验与观察: 控制粗轧的浪费
1. 用正态分布描述热轧机模型
2. 调整额定长度使浪费最小
.观察程序
zxy9_2.m
【 clf,clear, m=20;l=3;mu=3.5;sigma=0.3; mu1=3.5;sigma1=0.9;
a=0; b=max([mu+4*sigma,mu1+4*sigma1]); %设定坐标的范围。
subplot(2,2,1),
for k=1:m
rand('seed',1), x=normrnd(mu,sigma);
plot([0,x],[k,k],'linewidth',5),hold on, axis([0 mu+4*0.6 -2 m+5])
end
plot([l,l],[-2,m+5],'r-');hold on,axis([a b -2 m+5])
subplot(2,2,3),
t=linspace(a,b,50);f=normpdf(t,mu,sigma);
y0=max(f)+0.1;plot(t,f),hold on,axis([a b 0 y0]),
plot([l,l],[0,y0],'r-',[mu,mu],[0,y0],'r:');hold on
t=linspace(a,b,30);f=normpdf(t,mu,sigma);plot(t,f)
subplot(2,2,2)
for k=1:m
rand('seed',1),x=normrnd(mu1,sigma1);
plot([0,x],[k,k],'linewidth',5),hold on,axis([a b -2 m+5])
end
plot([l,l],[-2,m+5],'r-');hold on
subplot(2,2,4),
t=linspace(a,b,50);f=normpdf(t,mu1,sigma1);y0=max(f)+0.1;
plot(t,f),hold on,axis([a b 0 y0])
plot([l,l],[0,y0],'r-',[mu1,mu1],[0,y0],'r:');hold on 】
下面是观察与实验程序zxy9_3.m，对照上一节给出的算法，这一程序也是不难理解的。
zxy9_3.m
【 clf,clear
l=2;sigma=0.2;
n=10000;m=50;a=2.2;b=3;
mu=linspace(a,b,m);
for k=1:m
x=normrnd(mu(k),sigma,1,n); %模拟轧钢。
k_chenpin=find(x>=l); %求可轧制的成品材的下标。
k_feipin=find(x=(x(k,l)-h1)&x=(y(k,l)-h2)&y(j)<=(y(k,l)+h2));
z(k,l)=length(h)/n;
end
end
mu=[170 0.36*170]; %计算二维正态分布密度。
c=[4.5^2 0.36*4.5^2; 0.36*4.5^2, (0.36*4.5)^2+7^2];
detc=det(c);
for k=1:m
for l=1:m
u=[x(k,l),y(k,l)]'-mu';s=2*pi*sqrt(detc);
z1(k,l)=s^-1*exp((-0.5)*u'*c^-1*u);
end
end
subplot(2,2,1),surf(x,y,z),axis([a b c d 0 0.15]),
set(gca,'fontsize',15);
subplot(2,2,2),contour(x,y,z),axis([a b c d ]),axis('equal'),
subplot(2,2,3),surf(x,y,z1),axis([a b c d 0 0.005]),
subplot(2,2,4),contour(x,y,z1),axis([a b c d ]),axis('equal')】
3.用线性回归方法确定正常体重标准
观察：
◆（1）参考程序 zxy9_4.m中的模拟部分，。
【 alpha=0.01; polytool(x',y',1,alpha) 】
◆（3) 用 多元线性回归指令regress做体重预测。
◆对模拟的100对身高体重数据，运行下面的程序，了解指令regress 的用法。
【 [b,bint,r,rint,stats] = regress(y',[ones(100,1) x']);
b,bint,stats
rcoplot(r,rint),set(gca,'fontsize',15) 】
9.3实验（ⅲ）参数估计和假设检验
9.3.1实验与观察: 极大似然估计
1.极大似然估计原理: 如何决定废品率？
观察
◆（1）
【 clear,p=0.04; n=10; %设定产品的档次，抽样次数。
for k=1:n %抽样n次。
r(k)=rand(1,1); %产生随机数。
if r(k)p
x(k)=0; %抽样得到正品 。
end
end
i=[1:n];[i;x] 】
◆ （2）
2.实验观察的参考程序
实验观察的主要程序为zxy9_5.m，其源代码如下，该程序是不难读懂的。
zxy9_5.m
【 clear,clf,n=50;
p=0.04;
rand('seed',1),r=rand(1,n);
k=+find(r<=p);
x=zeros(1,n);
x(k)=ones(1,length(k));
p_estimate=sum(x)/n,
t=[0.01 0.02 0.03 0.04 0.05 0.06];
%t=linspace(0.01,0.5,40)
lt=t.^sum(x).*(1-t).^(n-sum(x));
%lnlt=sum(x)*log(t)+(n-sum(x))*log(1-t);
set(gca,'fontsize',16),
plot(t,lt,''),
[lmax,i]=max(lt);
tmax=t(i);
text(t(i),lmax,['\fontsize{16}\leftarrow\it lmax=' ,num2str(lmax)])
text(t(i),0.95*lmax,['\fontsize{16}\it{ at p}= ',num2str(tmax)])
xlabel('p','fontsize',16);ylabel('l(p)','fontsize',16); 】
9.3.2 应用、思考与练习
1. 用matlab符号演算求解极大似然估计
◆ 练习：用matlab符号演算指令求解9.3.1节中废品问题的似然方程获得极大似然估计。
【 syms p %未知参数为p,所以作为符号变量处理，用syms指令说明。
clear,clf,n=50; %产生50个样本。
p=0.04; %设定真实参数。
x=zeros(1,n); %令x全为0 。
rand('seed',1),r=rand(1,n);
k=+find(r<=p); %找出废品的下标。
x(k)=ones(1,length(k)); %在废品下标处改x为1；x为50个样本观察值。
%观察似然函数和似然方程的一般表达式。
l=sym('p^sx*(1-p)^(n-sx)') %正确写出似然函数，l是符号p的函数。
likely_equ=diff(l,'p') %对p进行符号求导，得到似然方程。
%观察含具体数值的似然函数和似然方程
sx=sum(x), p='p'; %代入sx的具体数值。
lp=subs(l) %将具体的数值代入似然函数中 。
likely_equ=diff(lp,'p') %求似然方程 。 】
【 %求似然方程的符号解，得到极大似然估计
s=solve('p^sx*sx/p*(1-p)^(n-sx)-p^sx*(1-p)^(n-sx)*(n-sx)/(1-p)=0','p')
sx,n %看看具体的数值。
sp=subs(s) %对已获得的样本，观察极大似然估计的具体数值。 】
2. 水库入库径流的分布估计
7月上旬径流数据
356 258 222 208 163 342 501 501 782 225 630 2305
931 485 503 501 422 101 280 1807 922 390 466 211
922 444 233 370 788 802 219 524 470 1097 1160 702
566 222 630
7月中旬径流数据
98 262 117 1687 291 1318 292 716 254 519 270 273
275 274 374 147 345 70 940 440 2839 141 699 324
900 311 870 596 187 2231 111 949 303 888 328 459
370 1360 1320
7月下旬径流数据
69 133 392 596 4518 1051 336 867 541 1733 149 266
324 1365 891 918 1751 219 513 438 1033 1217 1290 247 2360 1023 453 1622 1272 1383 1217 1530 1724 703 299 638
548 1200 1220
zxy9_6.m
【 clear,clf,
q=[ 98 262 117 1687 291 1318 292 716 254 519 270 273 275 274 374 147 345 70 940 440 2839 141 699 324 900 311 870 596 187 2231 111 949 303 888 328 459 370 1360 1320];
m=50;
as=linspace(0.6561, 2.2477,m);
bs=linspace(215.4687, 637.4421,m);
[a b]=meshgrid(as,bs);
for i=1:m
for j=1:m
ax=a(i,j);bx=b(i,j);
y(i,j) = sum(log(gampdf(q,ax,bx)));
end
end
[y0,s]=max(y); [y00,s0]=max(y0);
a_est=a(s0,s(s0));b_est=b(s0,s(s0)); meshc(a,b,y),hold on
plot3(a_est,b_est,y00,'.','markersize',25);
text(a_est,b_est,y00+0.06*abs(y00),['(a_e_s_t,b_e_s_t,l_m_a_x)=']);
text(a_est,b_est,y00+0.03*abs(y00),['(',num2str(a_est),',',num2str(b_est),',',num2str(y00),')'])
figure(2)
[h,n0]=hist(q,5); h0=max(h); a=min(q)-100;b=max(q);
x=linspace(a,b,30); hist(q,5); hold on
y = gampdf(x,a_est,b_est);
plot(x,h0*y/max(y),'r-','linewidth',5); 】
3. 数学建模竞赛: 零件的参数设计
zxy9_7.m
【 clear,clf, n=1000; a=0.01/3;b=0.05/3;c=0.1/3;
%x=[0.075,0.375,0.125,0.1198,1.252,12.5,0.77];%a=[b b b c c b b];
x=[0.1 0.3 0.1 0.1 1.5 16 0.75]; %设置标定值。
a=[b c c c c c b]; %设置容差。
x1=normrnd(x(1),a(1)*x(1),n,1); %模拟零件参数。
x2=normrnd(x(2),a(2)*x(2),n,1);x3=normrnd(x(3),a(3)*x(3),n,1);
x4=normrnd(x(4),a(4)*x(4),n,1);x5=normrnd(x(5),a(5)*x(5),n,1);
x6=normrnd(x(6),a(6)*x(6),n,1);x7=normrnd(x(7),a(7)*x(7),n,1);
y=z9_5fun(x1,x2,x3,x4,x5,x6,x7); %模拟y的样本。
k=find(imag(y)==0);y1=y(k); %如果产生了复数，剔除复数。
histfit(y1) %直方图和正态密度拟合。
dh=0.001; %求数值导数。
for i=1:7
for j=1:7
xx(:,j)=x(j)*ones(2,1);
end
xx(:,i)=linspace(x(i),x(i)+dh,2)';
f(:,i)=z9_5fun(xx(:,1),xx(:,2),xx(:,3),xx(:,4),xx(:,5),xx(:,6),xx(:,7));
g(:,i)=diff(f(:,i),1)/dh; %求数值导数。
end
mu=f(1,1),ey=mean(y1), %均值对比。
r=(a.*x).^2;
sigma=sqrt(sum((g.^2).*r)),dy=std(y1), %均方差值对比.
[h,p,ci]=ttest(y1,mu,0.01,0) %均值的t－检验 。 】
zxy9_6fun.m（计算函数的子程序）
【 function y=zxy9_5fun(x1,x2,x3,x4,x5,x6,x7)
y1=174.42*(x1./x5).*(x3./(x2-x1)).^0.85;
y2=(1-0.36*(x4./x2).^(-0.56)).^(1.5) ;
y3=(1-2.62*y2.*(x4./x2).^1.16)./(x6.*x7);
y=y1.*sqrt(y3);
200多个matlab经典教程和matlab论文请查看：matlab教程
到2025年，传统的计算技术将遭遇数字瓶颈
Dell发布2017首款新品：可翻转的XPS 13二合一笔电
MAX6033A电压基准解决方案
模拟电路和数字电路有何区别？
存储安全性必须在整体安全策略中得到应有的重视
matlab概率统计实验
瑞萨H3和骁龙665哪个好？
学生党蓝牙耳机推荐，性价比超高的运动蓝牙耳机
空间将成为最适合数据中心的冷却机
智能可穿戴设备——智能手表技术亮点解析
UHF RFID系统应用常见问题的认识与分析
SMT贴片加工为何如此受欢迎，它的特点是什么
让开发事半功倍！普华永道通过 Power Pages 实现效率升级
维修服务管理规范
异质结HIT电池的优缺点
S7-200以太网通讯处理器桥接型在纺机设备网络中的实际应用
双十一，论文分享
2019年比特币三大分支技术都有怎样的改变
小米12pro怎么样值得入手吗
荧光定量PCR仪的特点有哪些
9.1 实验（i）：galton钉板试验
9.1.1 实验与观察: galton钉板模型和二项分布
1. 动画模拟calton钉板试验
【 rand('seed',1), u=rand(1,6) 】
【 rand('seed',2),u=rand(1,6) 】
观察程序zxy9_1.m
【 clear，clf, m=100;n=5;y0=2; %设置参数。
ballnum=zeros(1,n+1);
p=0.5;q=1-p;
for i=n+1:-1:1 %创建钉子的坐标x,y 。
x(i,1)=0.5*(n-i+1);y(i,1)=(n-i+1)+y0;
for j=2:i
x(i,j)=x(i,1)+(j-1)*1;y(i,j)=y(i,1);
end
end
mm=moviein(m); % 动画开始，模拟小球下落路径。
for i=1:m
s=rand(1,n); %产生n个随机数。
xi=x(1,1);yi=y(1,1);k=1;l=1; % 小球遇到第一个钉子。
for j=1:n
plot(x(1:n,:),y(1:n,:),'o',x(n+1,:),y(n+1,:),'.-'), % 画钉子的位置 。
axis([-2 n+2 0 y0+n+1]),hold on
k=k+1; % 小球下落一格 。
if s(j)>p
l=l+0; %小球左移 。
else
l=l+1; %小球右移 。
end
xt=x(k,l);yt=y(k,l); %小球下落点的坐标。
h=plot([xi,xt],[yi,yt]);axis([-2 n+2 0 y0+n+1]) %画小球运动轨迹。
xi=xt;yi=yt;
end
ballnum(l)=ballnum(l)+1; %计数。
ballnum1=3*ballnum./m;
bar([0:n],ballnum1),axis([-2 n+2 0 y0+n+1]) %画各格子的频率 。
mm(i)=getframe; %存储动画数据。
hold off
end
movie(mm,1) %播放动画1次。 】
2. 用二项分布描述galton钉板模型
【 x = binornd(5,0.5,1,10) 】
3. 数学期望和平均收益
【 n=5;p=0.5; m=5000;
rand('seed',5); r = binornd(n,p,1,m); %模拟投球。
f=[5, 1, 0.2, 0.2, 1, 5]; %奖品的价值向量。
s=[];
for k=1:n+1 %计算第0～n格奖品价值。
u=[];
u=find(r==(k-1)); %计算落入第k-1格的小球下标，并存于向量u中。
s=[f(k)*length(u),s]; %计算相应的奖品价值，并存于向量s中。
end
mean_return=sum(s)/m %计算一次抽奖的平均回报 。 】
{ mean_return = 0.7506 }
【 f=[5, 1, 0.2, 0.2, 1, 5]; x=[0:1:n]; pu=binopdf(x,n,p);ef=sum(f.*pu) 】
9.1.2 应用、思考和练习
1. 二项分布的应用模型
◆电力供应问题 。
提示：
【 x=linspace(90,160,10); r = binornd(200,0.6,1,960); hist(r,x) ;
[n,x]=hist(r,x); 】
◆ 废品问题。
下面的程序可作为参考 。
【 clear,clf ,n=50;r=linspace(0,1,n);s=linspace(0,10,n);
[r,s]=meshgrid(r,s); x=[0 1 2 3 4 5];p = binopdf(x,5,0.5);
for k=1:n
f=[s(k,l) s(k,l)*r(k,l) s(k,l)*r(k,l)^2 s(k,l)*r(k,l)^2 s(k,l)*r(k,l) s(k,l)];
for l=1:n
z(k,l)=sum(f.*p);
if z(k,l)>=1 z(k,l)=nan; end
end
end
contour(r,s,z) 】
3. 单服务台定长服务时间排队系统的计算机模拟
4. 随机变量的模拟：逆概率法
图9.4 逆概率法模拟随机数原理。
【 clf
mu=0;sigma=1;a=mu-4*sigma;b=mu+4*sigma;
t=linspace(a,b,30);
f1=normcdf(t,mu,sigma);
for i=1:3
hold on
u=rand(1);
f=norminv(u,mu,sigma);
plot(t,f1,[0 0],[0,b]),hold on,
plot(0,u,'ro'),plot([0,f],[u,u]);plot([f,f],[u,0]);
plot(f,0,'ro'), axis([-4 4 0 1])
end 】
9.2 实验(ⅱ) ：热轧机的调整
9.2.1实验与观察: 控制粗轧的浪费
1. 用正态分布描述热轧机模型
2. 调整额定长度使浪费最小
.观察程序
zxy9_2.m
【 clf,clear, m=20;l=3;mu=3.5;sigma=0.3; mu1=3.5;sigma1=0.9;
a=0; b=max([mu+4*sigma,mu1+4*sigma1]); %设定坐标的范围。
subplot(2,2,1),
for k=1:m
rand('seed',1), x=normrnd(mu,sigma);
plot([0,x],[k,k],'linewidth',5),hold on, axis([0 mu+4*0.6 -2 m+5])
end
plot([l,l],[-2,m+5],'r-');hold on,axis([a b -2 m+5])
subplot(2,2,3),
t=linspace(a,b,50);f=normpdf(t,mu,sigma);
y0=max(f)+0.1;plot(t,f),hold on,axis([a b 0 y0]),
plot([l,l],[0,y0],'r-',[mu,mu],[0,y0],'r:');hold on
t=linspace(a,b,30);f=normpdf(t,mu,sigma);plot(t,f)
subplot(2,2,2)
for k=1:m
rand('seed',1),x=normrnd(mu1,sigma1);
plot([0,x],[k,k],'linewidth',5),hold on,axis([a b -2 m+5])
end
plot([l,l],[-2,m+5],'r-');hold on
subplot(2,2,4),
t=linspace(a,b,50);f=normpdf(t,mu1,sigma1);y0=max(f)+0.1;
plot(t,f),hold on,axis([a b 0 y0])
plot([l,l],[0,y0],'r-',[mu1,mu1],[0,y0],'r:');hold on 】
下面是观察与实验程序zxy9_3.m，对照上一节给出的算法，这一程序也是不难理解的。
zxy9_3.m
【 clf,clear
l=2;sigma=0.2;
n=10000;m=50;a=2.2;b=3;
mu=linspace(a,b,m);
for k=1:m
x=normrnd(mu(k),sigma,1,n); %模拟轧钢。
k_chenpin=find(x>=l); %求可轧制的成品材的下标。
k_feipin=find(x=(x(k,l)-h1)&x=(y(k,l)-h2)&y(j)<=(y(k,l)+h2));
z(k,l)=length(h)/n;
end
end
mu=[170 0.36*170]; %计算二维正态分布密度。
c=[4.5^2 0.36*4.5^2; 0.36*4.5^2, (0.36*4.5)^2+7^2];
detc=det(c);
for k=1:m
for l=1:m
u=[x(k,l),y(k,l)]'-mu';s=2*pi*sqrt(detc);
z1(k,l)=s^-1*exp((-0.5)*u'*c^-1*u);
end
end
subplot(2,2,1),surf(x,y,z),axis([a b c d 0 0.15]),
set(gca,'fontsize',15);
subplot(2,2,2),contour(x,y,z),axis([a b c d ]),axis('equal'),
subplot(2,2,3),surf(x,y,z1),axis([a b c d 0 0.005]),
subplot(2,2,4),contour(x,y,z1),axis([a b c d ]),axis('equal')】
3.用线性回归方法确定正常体重标准
观察：
◆（1）参考程序 zxy9_4.m中的模拟部分，。
【 alpha=0.01; polytool(x',y',1,alpha) 】
◆（3) 用 多元线性回归指令regress做体重预测。
◆对模拟的100对身高体重数据，运行下面的程序，了解指令regress 的用法。
【 [b,bint,r,rint,stats] = regress(y',[ones(100,1) x']);
b,bint,stats
rcoplot(r,rint),set(gca,'fontsize',15) 】
9.3实验（ⅲ）参数估计和假设检验
9.3.1实验与观察: 极大似然估计
1.极大似然估计原理: 如何决定废品率？
观察
◆（1）
【 clear,p=0.04; n=10; %设定产品的档次，抽样次数。
for k=1:n %抽样n次。
r(k)=rand(1,1); %产生随机数。
if r(k)p
x(k)=0; %抽样得到正品 。
end
end
i=[1:n];[i;x] 】
◆ （2）
2.实验观察的参考程序
实验观察的主要程序为zxy9_5.m，其源代码如下，该程序是不难读懂的。
zxy9_5.m
【 clear,clf,n=50;
p=0.04;
rand('seed',1),r=rand(1,n);
k=+find(r<=p);
x=zeros(1,n);
x(k)=ones(1,length(k));
p_estimate=sum(x)/n,
t=[0.01 0.02 0.03 0.04 0.05 0.06];
%t=linspace(0.01,0.5,40)
lt=t.^sum(x).*(1-t).^(n-sum(x));
%lnlt=sum(x)*log(t)+(n-sum(x))*log(1-t);
set(gca,'fontsize',16),
plot(t,lt,''),
[lmax,i]=max(lt);
tmax=t(i);
text(t(i),lmax,['\fontsize{16}\leftarrow\it lmax=' ,num2str(lmax)])
text(t(i),0.95*lmax,['\fontsize{16}\it{ at p}= ',num2str(tmax)])
xlabel('p','fontsize',16);ylabel('l(p)','fontsize',16); 】
9.3.2 应用、思考与练习
1. 用matlab符号演算求解极大似然估计
◆ 练习：用matlab符号演算指令求解9.3.1节中废品问题的似然方程获得极大似然估计。
【 syms p %未知参数为p,所以作为符号变量处理，用syms指令说明。
clear,clf,n=50; %产生50个样本。
p=0.04; %设定真实参数。
x=zeros(1,n); %令x全为0 。
rand('seed',1),r=rand(1,n);
k=+find(r<=p); %找出废品的下标。
x(k)=ones(1,length(k)); %在废品下标处改x为1；x为50个样本观察值。
%观察似然函数和似然方程的一般表达式。
l=sym('p^sx*(1-p)^(n-sx)') %正确写出似然函数，l是符号p的函数。
likely_equ=diff(l,'p') %对p进行符号求导，得到似然方程。
%观察含具体数值的似然函数和似然方程
sx=sum(x), p='p'; %代入sx的具体数值。
lp=subs(l) %将具体的数值代入似然函数中 。
likely_equ=diff(lp,'p') %求似然方程 。 】
【 %求似然方程的符号解，得到极大似然估计
s=solve('p^sx*sx/p*(1-p)^(n-sx)-p^sx*(1-p)^(n-sx)*(n-sx)/(1-p)=0','p')
sx,n %看看具体的数值。
sp=subs(s) %对已获得的样本，观察极大似然估计的具体数值。 】
2. 水库入库径流的分布估计
7月上旬径流数据
356 258 222 208 163 342 501 501 782 225 630 2305
931 485 503 501 422 101 280 1807 922 390 466 211
922 444 233 370 788 802 219 524 470 1097 1160 702
566 222 630
7月中旬径流数据
98 262 117 1687 291 1318 292 716 254 519 270 273
275 274 374 147 345 70 940 440 2839 141 699 324
900 311 870 596 187 2231 111 949 303 888 328 459
370 1360 1320
7月下旬径流数据
69 133 392 596 4518 1051 336 867 541 1733 149 266
324 1365 891 918 1751 219 513 438 1033 1217 1290 247 2360 1023 453 1622 1272 1383 1217 1530 1724 703 299 638
548 1200 1220
zxy9_6.m
【 clear,clf,
q=[ 98 262 117 1687 291 1318 292 716 254 519 270 273 275 274 374 147 345 70 940 440 2839 141 699 324 900 311 870 596 187 2231 111 949 303 888 328 459 370 1360 1320];
m=50;
as=linspace(0.6561, 2.2477,m);
bs=linspace(215.4687, 637.4421,m);
[a b]=meshgrid(as,bs);
for i=1:m
for j=1:m
ax=a(i,j);bx=b(i,j);
y(i,j) = sum(log(gampdf(q,ax,bx)));
end
end
[y0,s]=max(y); [y00,s0]=max(y0);
a_est=a(s0,s(s0));b_est=b(s0,s(s0)); meshc(a,b,y),hold on
plot3(a_est,b_est,y00,'.','markersize',25);
text(a_est,b_est,y00+0.06*abs(y00),['(a_e_s_t,b_e_s_t,l_m_a_x)=']);
text(a_est,b_est,y00+0.03*abs(y00),['(',num2str(a_est),',',num2str(b_est),',',num2str(y00),')'])
figure(2)
[h,n0]=hist(q,5); h0=max(h); a=min(q)-100;b=max(q);
x=linspace(a,b,30); hist(q,5); hold on
y = gampdf(x,a_est,b_est);
plot(x,h0*y/max(y),'r-','linewidth',5); 】
3. 数学建模竞赛: 零件的参数设计
zxy9_7.m
【 clear,clf, n=1000; a=0.01/3;b=0.05/3;c=0.1/3;
%x=[0.075,0.375,0.125,0.1198,1.252,12.5,0.77];%a=[b b b c c b b];
x=[0.1 0.3 0.1 0.1 1.5 16 0.75]; %设置标定值。
a=[b c c c c c b]; %设置容差。
x1=normrnd(x(1),a(1)*x(1),n,1); %模拟零件参数。
x2=normrnd(x(2),a(2)*x(2),n,1);x3=normrnd(x(3),a(3)*x(3),n,1);
x4=normrnd(x(4),a(4)*x(4),n,1);x5=normrnd(x(5),a(5)*x(5),n,1);
x6=normrnd(x(6),a(6)*x(6),n,1);x7=normrnd(x(7),a(7)*x(7),n,1);
y=z9_5fun(x1,x2,x3,x4,x5,x6,x7); %模拟y的样本。
k=find(imag(y)==0);y1=y(k); %如果产生了复数，剔除复数。
histfit(y1) %直方图和正态密度拟合。
dh=0.001; %求数值导数。
for i=1:7
for j=1:7
xx(:,j)=x(j)*ones(2,1);
end
xx(:,i)=linspace(x(i),x(i)+dh,2)';
f(:,i)=z9_5fun(xx(:,1),xx(:,2),xx(:,3),xx(:,4),xx(:,5),xx(:,6),xx(:,7));
g(:,i)=diff(f(:,i),1)/dh; %求数值导数。
end
mu=f(1,1),ey=mean(y1), %均值对比。
r=(a.*x).^2;
sigma=sqrt(sum((g.^2).*r)),dy=std(y1), %均方差值对比.
[h,p,ci]=ttest(y1,mu,0.01,0) %均值的t－检验 。 】
zxy9_6fun.m（计算函数的子程序）
【 function y=zxy9_5fun(x1,x2,x3,x4,x5,x6,x7)
y1=174.42*(x1./x5).*(x3./(x2-x1)).^0.85;
y2=(1-0.36*(x4./x2).^(-0.56)).^(1.5) ;
y3=(1-2.62*y2.*(x4./x2).^1.16)./(x6.*x7);
y=y1.*sqrt(y3);
200多个matlab经典教程和matlab论文请查看：matlab教程
到2025年，传统的计算技术将遭遇数字瓶颈
Dell发布2017首款新品：可翻转的XPS 13二合一笔电
MAX6033A电压基准解决方案
模拟电路和数字电路有何区别？
存储安全性必须在整体安全策略中得到应有的重视
matlab概率统计实验
瑞萨H3和骁龙665哪个好？
学生党蓝牙耳机推荐，性价比超高的运动蓝牙耳机
空间将成为最适合数据中心的冷却机
智能可穿戴设备——智能手表技术亮点解析
UHF RFID系统应用常见问题的认识与分析
SMT贴片加工为何如此受欢迎，它的特点是什么
让开发事半功倍！普华永道通过 Power Pages 实现效率升级
维修服务管理规范
异质结HIT电池的优缺点
S7-200以太网通讯处理器桥接型在纺机设备网络中的实际应用
双十一，论文分享
2019年比特币三大分支技术都有怎样的改变
小米12pro怎么样值得入手吗
荧光定量PCR仪的特点有哪些