第12次实验作业

实验报告（二）
院（系）理学院课程名称：数学模型日期2012-05-30 班级A0921 学号01 实验室文理楼209 专业信息与计算科学姓名计算机号12
实验
名称
概率论实验成绩评定
所用
软件
matlab 教师签名实
验
目的或要求1.了解和掌握使用蒙特卡罗方法
2.掌握回归分析
实验步骤、心得体会实验题目
1.用蒙特卡罗方法计算下面两个积分
(1)22
d
x x
⎰(2)10sin d.
x
x
x
⎰
解：（1）输入
clear
n=1000000;m=20;
b=ones(1);
for k=1:m
a=rand(n,2);
x(:,1)=2*a(:,1);
x(:,2)=4*a(:,2);
count=0;
for i=1:n
if x(i,1)^2>=x(i,2); count=count+1;
end
end
b(k)=count/n;
end
disp(8*mean(b))
输出
2.6659 即
2
20
2.6659x dx ≈⎰
（2）输入
clear
n=1000000; x=rand(n,2); count=0;
for i=1:n
if (sin(x(i,1)))/x(i,1)>=x(i,2); count=count+1; end
end
disp(count/n); 输出
0.9462 即9462.0sin 1
0≈⎰
dx x
x
2.研究表明，合金的强度y 与含碳量x 存在某种关系。

现有一批数据，试研究着两者之间的关系。

x 0.10 0.11 0.12 0.13 0.14 0.15 y 41.0 42.5 45.0 45.5 45.0 47.5 x 0.16 0.17 0.18 0.20 0.22 0.24 y
49.0
51.0
50.0
55.5
57.5
59.5
解：作出散点图
x=[0.10 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.20 0.22 0.24];
y=[41.0 42.5 45.0 45.5 45.0 47.5 49.0 51.0 50.0 55.5 57.5 59.5];
A=polyfit(x,y,1) z=polyval(A,x);
plot(x,y,'r+',x,z,'b') grid on
根据散点图我们相信合金的强度y与含碳量x之间的关系是线性关系，下面我们进行一元回归分析，
m=length(x);
X=[ones(m,1),x'];
[b bint r rint stats]=regress(y',X,0.05)
figure(1),rcoplot(r,rint)
y1=b(1)+b(2)*x;
figure(2),plot(x,y,'r*',x,y1)
输出
b =
27.8381
132.7830
bint =
25.4860 30.1901
118.5649 147.0011
stats =
0.9774 432.9980 0.0000 0.8632
残差图如下：
在此残差图中，第9个数据点异常，将其踢除后重新检验，得残差图此时，第5个数据点异常，将其剔除，再次检验，得残差图
此时相应的数值为：
b =
28.1059
133.0882
bint =
26.5859 29.6259
123.9454 142.2311
stats =
1.0e+003 *
0.0010 1.1268 0.0000 0.0003
残差的置信区间都包含零点，说明回归模型较为理想。

所以，合金的强度y与含碳量x之间的关系是
=+
y x
133.088228.1059
关系图如下：
3.某地区车祸次数y（千次）与汽车拥有量x（万辆）的11年统计数据如下表
年度汽车量车祸次数年度汽车量车祸次数
1 35
2 166 7 529 227
2 37
3 153 8 577 238
3 411 177 9 641 268
4 441 201 10 692 268
5 462 21
6 11 743 274
6 490 208
（1）作出y和x的散点图
（2）如果从散点图大致看出y对x是线性的，试求线性回归方程
解：（1）程序
y=[166 153 177 201 216 208 227 238 268 268 274];
x=[352 373 411 441 462 490 529 577 641 692 743];
plot(x,y,'r+')
grid on
得到y和x的散点图如下
（2）根据散点图，我们相信y与x之间的关系是线性的，下面我们求线性回归方程y=[166 153 177 201 216 208 227 238 268 268 274];
x=[352 373 411 441 462 490 529 577 641 692 743];
m=length(x);
X=[ones(m,1),x'];
[b bint r rint stats]=regress(y',X,0.05)
figure(1),rcoplot(r,rint)
y1=b(1)+b(2)*x;
figure(2),plot(x,y,'r*',x,y1)
得到残差图如下
由此可知，第2个数据点异常，将其剔除，重新检验，得残差图
此时相应的值为
b =
69.0716
0.2908
bint =
39.2370 98.9062
0.2363 0.3453
stats =
0.9498 151.3285 0.0000 81.5805
残差的置信区间都包含零点，说明回归模型较为理想。

所以，线性回归方程为0.290869.0716
=+
y x
关系图如下
4.铅墙的高度远大于左右厚度。

设中子是垂直由左端进入铅墙，在铅墙中运行一个单位距离然后与一个铅原子碰撞。

碰撞后，任意改变方向，并继续运行一个单位后与另一个铅原子碰撞。

这样下去，如果中子在铅墙里消耗掉所有的能量或者从左端逸出就被视为中子被铅墙挡住，如果中子穿过铅墙由右端逸出就视为中子逸出。

如果铅墙厚度为5个单位，中子运行7个单位后能
量耗尽，求中子逸出几率。

解：建模
我们想要了解的是一次碰撞后，中子在x轴上行进的距离是多少，所以行进方向是正负θ的结果是一样的,我们就只考虑θ是正的情形. 由于中子运行的方向θ是随机的,我们用计算机抽取在0到π间均衡分布的随机数,模拟1000000个中子在铅墙里行进的情形,看看这些中子与铅原子碰撞7次后,有多少超过了铅墙的右端.
程序
clear,clc
n=1000000;m=0;
for i=1:n
x=1;t=2;
for k=1:7
ang=pi*rand;
x=x+cos(ang);
if x<0
t=0;
break
end
if x>5
t=1;
break
end
end
if t==1
m=m+1;
end
end
disp(m/n)
输出
0.0187
得出结果，中子逸出几率为0.0187
备注：本实验报告用于各学科与计算机应用相关课程的实验，务必按时完成。

不交此报告者，本次实验为“不合格”。

九江学院。