《概率统计的数值实验——t整理.ppt》由会员分享，可在线阅读，更多相关《概率统计的数值实验——t整理.ppt（92页珍藏版）》请在金锄头文库上搜索。
1、概率统计的数值实验 MATLAB在概率统计教学中的应用 崔明涛2012年10月11日 引言 而MATLAB软件具有简单易学 易操作和绘图功能强等特点 利用MATLAB软件的图形可视功能将概率统计的内容用图形表示出来 通过图形让学生加深理解 以达到事半功倍的效果 概率论与数理统计知识比较抽象 逻辑性较强 因此 建议让学生结合理论和公式推导 进行数值试验和相关调查 直观地感受数学概念和理论 从而提高学生解决实际问题的信心和能力 概率论 1 rand m n 生成m n的随机矩阵 每个元素都在 0 1 间 生成方式为均匀分布 2 randn m n 生成m n的随机矩阵 每个元素都在 0 1 间 生
2、成方式为正态分布 3 randperm m 生成一个1 m的随机整数排列 4 perms 1 n 生成一个1 n的全排列 共n 个 5 取整函数系列 1 fix x 截尾法取整 2 floor x 退一法取整 不超过x的最大整数 3 ceil x 进一法取整 floor x 1 4 round x 四舍五入法取整 6 unique a 合并a中相同的项 7 prod x 向量x的所有分量元素的积 一 MATLAB常用的与随机数产生相关的函数 示例 rand 1 生成一个 0 1 间的随机数ans 0 8147 rand 2 2 生成一个2 2阶 0 1 间的随机数矩阵ans 0 91340 0
3、9750 63240 2785 randperm 5 生成一个1 5的随机整数排列ans 41523 a 1242332 unique a ans 1234 例1随机投掷均匀硬币 观察国徽朝上与国徽朝下的频率 解 n 3000 100000000 m 0 fori 1 nt randperm 2 生成一个1 2的随机整数排列x t 1 生成一个0 1的随机整数排列y x 1 ify 0 m m 1 endendp1 m np2 1 p1 可见当时 解记事件为第i个人拿到自已枪 事件为第i个人没拿到自己枪 易知 又记为没有一个人拿到自己枪的概率 有乘法公式可知 例2某班有n个人 每人各有一支枪
4、这些枪外形一样 某次夜间紧急集合 若每人随机地取走一支枪 问没有一个人拿到自己枪的概率是多少 于是所以特别地 当n较大时 因此 可随机模拟出没有人拿到自己枪的频率 根据频率的稳定性 近似当做概率 然后去估计自然常数e 算法如下 1 产生n个随机数的随机序列 2 检验随机列与自然列是否至少有一个配对 3 对没有一个配对的序列进行累积p 4 重复1 2 3步m次 5 估计 具体程序及相关结果为 注 自然常数e 2 7183 m 40000 n 50 p 0 forj 1 mk 0 sui randperm n fori 1 nifsui i ik k 1 elsek k endendifk 0p
5、p 1 elsep p endende m pe 2 7313 设针与平行线的夹角为 针的中心与最近直线的距离为 针与平行线相交的充要条件是 则所求概率为故可得的近似计算公式 其中n为随机试验次数 m为针与平行线相交的次数 例3Buffon投针实验在画有许多间距为的等距平行线的白纸上 随机投掷一根长为的均匀直针 求针与平行线相交的概率 并计算的近似值 解 clear clfn 10000000 l 0 5 m 0 d 1 fori 1 nx l 2 sin rand 1 pi y rand 1 d 2 ifx ym m 1 endendp1 m npai 2 n l m d 例4在100个人的
6、团体中 不考虑年龄差异 研究是否有两个以上的人生日相同 假设每人的生日在一年365天中的任意一天是等可能的 那么随机找n个人 不超过365人 1 求这n个人生日各不相同的概率是多少 从而求这n个人中至少有两个人生日相同这一随机事件发生的概率是多少 2 近似计算在30名学生的一个班中至少有两个人生日相同的概率是多少 解 1 clear clfforn 1 100p0 n prod 365 1 365 n 1 365 n p1 n 1 p0 n endp1 ones 1 100 p0 n 1 100 plot n p0 n p1 xlabel 人数 ylabel 概率 legend 生日各不相同的
7、概率 至少两人生日相同的概率 axis 0100 0 11 199 gridon p1 30 0 7063 p1 60 0 9941 分析 在30名学生中至少两人生日相同的概率为70 63 下面进行计算机仿真 随机产生30个正整数 代表一个班30名学生的生日 然后观察是否有两人以上生日相同 当30个人中有两人生日相同时 输出 1 否则输出 0 如此重复观察100次 计算出这一事件发生的频率 2 clear clfn 0 form 1 100 做100次随机试验y 0 x 1 fix 365 rand 1 30 产生30个随机数fori 1 29 用二重循环寻找30个随机数中是否有相同数forj
8、 i 1 30ifx i x j y 1 break endendendn n y 累计有两人生日相同的试验次数endf n m 计算频率 f 0 6900f 0 7900f 0 6700f 0 7300f 0 7500f 0 6900f 0 7200f 0 6700f 0 6800 重复观察 数据如下 例5Galton钉板模型和二项分布Galton钉板试验是由英国生物统计学家和人类学家Galton设计的 故而得名 通过模拟Calton钉板试验 观察和体会二项分布概率分布列的意义 形象地理解DeMoivre Laplace中心极限定理 共15层小钉 高尔顿钉板试验 小球最后落入的格数 记小球向
9、右落下的次数为则 记小球向左落下的次数为则 符号函数 大于0返回1 小于0返回 1 等于0返回0 高尔顿 FrancisGalton 1822 1911 英国人类学家和气象学家 O 记 则 近似 高尔顿钉板试验 共15层小钉 模拟Galton钉板试验的步骤 1 确定钉子的位置 将钉子的横 纵坐标存储在两个矩阵X和Y中 2 在Galton钉板试验中 小球每碰到钉子下落时都具有两种可能性 设向右的概率为p 向左的概率为q 1 p 这里p 0 5 表示向左向右的机会是相同的 模拟过程如下 首先产生一均匀随机数u 这只需调用随机数发生器指令rand m n rand m n 指令 用来产生m n个 0
10、 1 区间中的随机数 并将这些随机数存于一个m n矩阵中 每次调用rand m n 的结果都会不同 如果想保持结果一致 可与rand seed s 配合使用 这里s是一个正整数 例如 rand seed 1 u rand 1 6 u 0 51290 46050 35040 09500 43370 7092而且再次运行该指令时结果保持不变 除非重设种子seed的值 如 rand seed 2 u rand 1 6 u 0 02580 92100 70080 19010 86730 4185这样结果才会产生变化 将 0 1 区间分成两段 区间 0 p 和 p 1 如果随机数u属于 0 p 让小球向
11、右落下 若u属于 p 1 让小球向左落下 将这一过程重复n次 并用直线连接小球落下时所经过的点 这样就模拟了小球从顶端随机地落人某一格子的过程 3 模拟小球堆积的形状 输入扔球次数m 例如m 50 100 500等等 计算落在第i个格子的小球数在总球数m中所占的比例 这样当模拟结束时 就得到了频率用频率反映小球的堆积形状 4 用如下动画指令制作动画 movien n 创建动画矩阵 制作动画矩阵数据 Getframe 拷贝动画矩阵 movie Mat m 播放动画矩阵m次 M文件如下 解 clear clf m 100 n 5 y0 2 设置参数ballnum zeros 1 n 1 p 0 5
12、 q 1 p fori n 1 1 1 创建钉子的坐标x yx i 1 0 5 n i 1 y i 1 n i 1 y0 forj 2 ix i j x i 1 j 1 1 y i j y i 1 endendmm moviein m 动画开始 模拟小球下落路径fori 1 ms rand 1 n 产生n个随机数xi x 1 1 yi y 1 1 k 1 l 1 小球遇到第一个钉子forj 1 nplot x 1 n y 1 n o x n 1 y n 1 画钉子的位置axis 2n 20y0 n 1 holdon k k 1 小球下落一格ifs j pl l 0 小球左移elsel l 1
13、小球右移endxt x k l yt y k l 小球下落点的坐标h plot xi xt yi yt axis 2n 20y0 n 1 画小球运动轨迹xi xt yi yt endballnum l ballnum l 1 计数ballnum1 3 ballnum m bar 0 n ballnum1 axis 2n 20y0 n 1 画各格子的频率mm i getframe 存储动画数据holdoffendmovie mm 1 播放动画一次 概率密度函数 pdf 求随机变量X在x点处的概率密度值 累积分布函数 cdf 求随机变量X在x点处的分布函数值 逆累积分布函数 inv 求随机变量X在
14、概率点处的分布函数反函数值 均值与方差计算函数 stat 求给定分布的随机变量X的数学期望E X 和方差var X 随机数生成函数 rnd 模拟生成指定分布的样本数据 二 MATLAB为常见自然概率分布提供了下列5类函数 具体函数的命名规则是 函数名 分布类型名称 函数类型名称 pdf cdf inv stat rnd 其中 分布类型名称如下 分布类型MATLAB名称 正态分布norm指数分布exp均匀分布unif 分布beta 分布gam对数正态分布lognrayleigh分布raylweibull分布weib二项分布binoPoisson分布poiss几何分布geo超几何分布hyge离散均
15、匀分布unid负二项分布nbin 例如 normpdf normcdf norminv normstat和normrnd分别是正态分布的概率密度 累积分布 逆累积分布 数字特征和随机数生成函数 关于这5类函数的语法 请详见有关书籍 快捷的学习可借助MATLAB的系统帮助 通过指令doc获得具体函数的详细信息 语法是doc 例6到某服务机构办事总是要排队等待的 设等待时间T是服从指数分布的随机变量 单位 分钟 概率密度为设某人一个月内要到此办事10次 若等待时间超过15分钟 他就离去 求 1 恰好有两次离去的概率 2 最多有两次离去的概率 3 至少有两次离去的概率 4 离去的次数占多数的概率 解
16、首先求任一次离去的概率 依题意设10次中离去的次数为X 则 p 1 expcdf 15 10 任一次离去的概率p1 binopdf 2 10 p 恰有两次离去的概率q binopdf 0 2 10 p p2 sum q 最多有两次离去的概率q binopdf 0 1 10 p p3 1 sum q 最少有两次离去的概率q binopdf 0 5 10 p p4 1 sum q 离去的次数占多数的概率 p 0 2231p1 0 2972p2 0 6073p3 0 6899p4 0 0112 例7某一急救中心在长度为t的时间间隔内收到的紧急呼救次数服从参数为t 2的泊松分布 而与时间间隔的起点无关 时间以小时计 求 1 在某一天中午12时至下午3时没有收到紧急呼救的概率 2 某一天中午12时至下午5时至少收到1次紧急呼救的概率 1 P1 poisscdf 0 3 2 P1 0 2231或者 P1 poisspdf 0 3 2 P1 0 2231中午12时到下午3时没有收到紧急呼救的概率为0 2231 2 P2 1 poisscdf 0 5 2 P2 0 9179中午12时至下午5时至少收到
