蒙特卡罗法(Monte Carlo ).doc

实验十二计算机仿真实验实验目的：1. 掌握全概率公式与贝叶斯公式；2. 了解计算机仿真方法；3. 了解蒙特卡罗法(Monte Carlo method), 具有初级编程能力.实验原理：全概率公式: 设A 1, A2, , An为两两互斥事件，B是A 1 + A2 + + An的子事件，则P(B)=P(A1)P(B|A1) + P(A2)P(B|A2) + + P(An)P(B|An).贝叶斯公式: P(Ak|B)= P(Ak )P(B| Ak )/P(B).计算机仿真: 就是在计算机上模拟各种实际系统的运行过程. 计算机仿真通常用来产生规定分布的随机变量.对于任意随机变量x ，其分布函数为F(x)，设 h = F(x )的分布函数为G(y)，则G(y) = Ph y= PF(x ) y= Px F -1( y)= y，这说明h 服从0，1的均匀分布.一般的编程语言都提供了均匀分布随机数发生器.应用随机数模拟试验的方法通常称为蒙特卡罗法(Monte Carlo method). 蒙特卡罗法不仅适用于处理随机性问题, 如存贮、排队、质量检验、市场营销、社会救急、生态竞争和传染病等问题；也可处理定性问题, 如计算多重积分、解积分方程及微分方程、解整数规划(特别是非线性整数规划)等.应用蒙特卡罗法解规划问题的基本思想是：先估计各个变量的大致取值范围，每次试验从中随机取出一个样本点，然后判断它是否为可行点. 若是则将其目标函数值与上一次的目标函数值相比较，记录下较优目标函数值与其样本点；否则重新抽样。直到试验次数达到指定值或可行点数达到指定值为止.实验内容：1. 设有两个口袋，甲袋中盛有两个白球，一个黑球，乙袋中盛有一个白球，两个黑球.由甲袋任取一个球放入乙袋，再从乙袋中取出一个球.若从乙袋中取出的球是白球，那么从甲袋中取出放入乙袋的球是白球还是黑球? 用计算机模拟上述过程1000次，问理论判断是正确的有多少次?2. 用计算机模拟随机变量x j (x) =的取值200次.3. 用计算机模拟随机变量x N (120,102 )的取值800次，并画出统计直方图.4. 应用蒙特卡罗法解非线性规划问题:max z = - 2x2 - y2 + xy + 8x + 3ys.t. 3 x + y = 10x0, y0.C语言简介(仅介绍本实验所用到的)1.标识符标识符是由程序员定义的单词，如函数名、变量名等. 标识符是由大小写字母、数字和下划线组成的，并以字母和下划线开始.2.关键字void (无值型) char (字符型) int (整数型) long (长整数型) float (浮点型) double (双精度浮点型)if (如果) else (或者)for (循环) while (循环) break (满足一定的条件终止循环)return (返回函数值)3.函数形式类型 函数名(参数)4.库函数简介double sin(double x) double cos(double x) double exp(double x) double log(double x)double sqrt(double x) double pow(double x, double y)double fabs(double x) int abs(int x)int printf(const char *format, .); 屏幕格式化输出函数FILE *fopen(const char *filename, const char *mode); 文件打开int fprintf(FILE *fp, const char *format, .); 文件格式化输出函数int fgetc(FILE *fp); 从文件中读出一字符int fclose(FILE *fp); 文件关闭void far setcolor(int color); 设置输出颜色void far line(int x1, int y1, int x2, int y2); 画直线void far rectangle(int left, int top, int right, int bottom); 画矩形unsigned far getpixel(int x, int y); 读出点(x, y)的颜色void far putpixel(int x, int y, int pixelcolor); 画点int random(int Num), 均匀产生0到Num-1中的一个随机数5.示例计算9！#includevoid main(void) int i; long n; n=1; for(i=1;i=9;i+)n*=i; printf(n9!=%ldn,n); getch();6. Turbo C(2.0) 编辑命令F3 录入文件F9 编译Ctrl+F9 运行Ctrl+KB 定义块首Ctrl+KK 定义块尾Ctrl+KC 块粘贴Ctrl+KV 块移动Ctrl+Y 删除当前行7.部分源程序程序LAB1_1.C 求出方程sinx - x = 1在( -2，2)内的近似根#include#includefloat f(float x) return sin(x)-x-1; void main(void) float r=0.618,x0=-2,x1=2,x; int n=0; while(1) n+; x=(1-r)*x0+r*x1; if(f(x)*f(x0)0.0)x1=x; else if(f(x)*f(x1)0.0)x0=x; else break; if(fabs(x1-x0)0.001)break; printf(n=%d, x=%fn,n,x); getch();程序LAB2_3.C 给出正态分布函数表#include#includefloat f(float x) return exp(-x*x/2); void main(void) float x,x0,x1=0.0,F=0,h=0.0001; long n=0; FILE *fp; int i=0,p; fp=fopen(x.c,w); for(x1=0.0;x10.04;x1+=0.001) n=0; F=0; x0=-10; for(x=x0+h;xx1;x+=h) if(n%2)F+=2*f(x); else F+=4*f(x); n+; F+=f(x0)+f(x1); F/=3; F*=h; F*=0.39894; printf(x=%5.3f, F=%6.4fn,x1,F); i+; p=F*10000; fprintf(fp,%d,p); if(i=10) i=0; fprintf(fp,n); if(p=9999)break; fclose(fp); getch();程序LAB3_2_1.C 解下列微分方程y- ytanx = secx，y(0) = 0, 并画出其图形:#include #include #include float f(float x,float y) return y*sin(x)/cos(x)+1/cos(x); void main(void) int i=DETECT,j; float x=0.0,y=0.0,h=0.005; char *str=0.00; initgraph(&i,&j, ); setviewport(0,0,639,479,1); cleardevice(); setbkcolor(BLUE); setcolor(WHITE); line(20,200,620,200); for(i=0;i10;i+) line(20+i*60,195,20+i*60,200); str0=48+3*i/10; str2=48+3*i%10; outtextxy(20+i*60,205,str); for(i=0;i600;i+) y=y+f(x,y)*h; x+=h; putpixel(i+20,200-y*10,GREEN); getch(); closegraph();程序LAB3_2_2.C 解下列微分方程组 x (0) = 100, y (0) = 8, 并画出其图形:#include #include #include float f1(float x,float y) return 2*x-0.25*x*y; float f2(float x,float y) return -y+0.01*x*y; void main(void) int i=DETECT,j; float x=99.0,y=7.9,h=0.015; initgraph(&i,&j, ); setviewport(0,0,639,479,1); cleardevice(); for(i=0;i600;i+) x=x+f1(x,y)*h; y=y+f2(x,y)*h; putpixel(i+20,(102.5-x)*60,GREEN); putpixel(i+20,(8.6-y)*600,WHITE); getch(); closegraph();程序LAB4.C 将矩阵化为行阶梯型,化为行最简型#include#include#define MAXR 20#define MAXC 40/*解线性方程组,以下是其增广矩阵*/float M_BMAXRMAXC=1,-2,2,1,-3,2,1,1,-2,-1,3,4,0,-5,1;void f1(int m,int n) int i,j,r=0,c=0; float x0; printf(n以下将矩阵化为行阶梯型n); for(i=0;im;i+) for(j=0;jn;j+)printf(%8.2f,M_Bij); printf(n); printf(按任一健继续.n); getch(); while(cn) for(i=r;i=0.0001)break; if(im) if(i!=r)for(j=0;jn;j+)x0=M_Bij;M_Bij=M_Brj;M_Brj=x0; x0=M_Brc; for(j=0;jn;j+)M_Brj/=x0; for(i=r+1;im;i+) x0=M_Bic;for(j=0;jn;j+)M_Bij-=x0*M_Brj; r+; for(i=0;im;i+) for(j=0;jn;j+)printf(%8.2f,M_Bij); printf(n); printf(按任一健继续.n); getch(); c+; if(r=m)break; printf(n以下将行阶梯型化为行最简型n); while(r) r-; for(j=0;j=0.0001)break; c=j; for(i=0;ir;i+) x0=M_Bic; for(j=0;jn;j+)M_Bij-=x0*M_Brj; for(i=0;im;i+) for(j=0;jn;j+)printf(%8.2f,M_Bij); printf(n); printf(按任一健继续.n); getch(); printf(完毕,按任一健退出.n); getch();void f2(int n) int i,j; for(i=0;in;i+)for(j=0;jn;j+) if(i!=j)M_Bin+j=0; else M_Bin+j=1; f1(n,2*n);void main(void) f1(3,5); 程序LAB6_3.C 模拟产生服从N(120,400)分布的随机变量800次, 并画出统计直方图#include #include #include #include int Np4000=/*下列数据为标准正态分布函数值10000, 由程序LAB2_3.C 产生*/5000,5003,5007,5011,5015,5019,5023,5027,5031,5035,;int f(int p) int i; for(i=0;i=p)break; return i;void main(void) int i,j,x; float F,X800; randomize(); printf(nn); for(i=0;i800;i+) x=random(10000); if(x5000)F=-f(10000-x); else F=f(x); F*=0.001; Xi=120+20*F; printf(%8.2f,Xi); getch(); i=DETECT; initgraph(&i,&j, ); setviewport(0,0,639,479,1); cleardevice(); setbkcolor(BLUE); setcolor(WHITE); line(20,400,620,400); j=20; for(F=90;F150;F+=3) x=0; for(i=0;i=F&XiF+5)x+=4; rectangle(j,400-x,j+30,400); j+=30; getch(); closegraph();程序2000A01.C 统计文件2000A1.txt中ATCG的个数