蒙特卡罗法(Monte Carlo method)

计算机仿真实验

实验目的：

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).

计算机仿真: 就是在计算机上模拟各种实际系统的运行过程. 计算机仿真通常用来产生规定分布的随机变量.

对于任意随机变量  ，其分布函数为F(x)，设   = F(   )的分布函数为G(y)，则G(y) = P{  ≤y}= P{F(   ) ≤y}= P{   ≤F  1( y)}= y，这说明  服从[0，1]的均匀分布.一般的编程语言都提供了均匀分布随机数发生器.

应用随机数模拟试验的方法通常称为蒙特卡罗法(Monte Carlo method). 蒙特卡罗法不仅适用于处理随机性问题, 如存贮、排队、质量检验、市场营销、社会救急、生态竞争和传染病等问题；也可处理定性问题, 如计算多重积分、解积分方程及微分方程、解整数规划(特别是非线性整数规划)等.

应用蒙特卡罗法解规划问题的基本思想是：先估计各个变量的大致取值范围，每次试验从中随机取出一个样本点，然后判断它是否为可行点. 若是则将其目标函数值与上一次的目标函数值相比较，记录下较优目标函数值与其样本点；否则重新抽样。直到试验次数达到指定值或可行点数达到指定值为止.

实验内容：

1. 设有两个口袋，甲袋中盛有两个白球，一个黑球，乙袋中盛有一个白球，两个黑球.由甲袋任取一个球放入乙袋，再从乙袋中取出一个球.若从乙袋中取出的球是白球，那么从甲袋中取出放入乙袋的球是白球还是黑球? 用计算机模拟上述过程1000次，问理论判断是正确的有多少次?

2. 用计算机模拟随机变量   ～  (x) =的取值200次.

3. 用计算机模拟随机变量  ～N (120,102 )的取值800次，并画出统计直方图.

4. 应用蒙特卡罗法解非线性规划问题:

max  z =   2x2   y2 + xy + 8x + 3y

s.t.  3 x + y = 10

x≥0, y≥0.

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！

#include

void main(void)

{

 int i;

 long n;

 n=1;

for(i=1;i<=9;i++)n*=i;

 printf("\\n9!=%ld\\n",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

#include

float 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=%f\\n",n,x);

  getch();

}

程序LAB2_3.C 给出正态分布函数表

#include

#include

float 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

for(x1=0.0;x1<0.04;x1+=0.001){

    n=0;    F=0;    x0=-10;

for(x=x0+h;x      if(n%2)F+=2*f(x);

      else F+=4*f(x);

      n++;

    }

    F+=f(x0)+f(x1);

    F/=3;    F*=h;    F*=0.394;

    printf("x=%5.3f,  F=%6.4f\\n",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;i<10;i++){

    line(20+i*60,195,20+i*60,200);

    str[0]=48+3*i/10;

    str[2]=48+3*i%10;

    outtextxy(20+i*60,205,str);

  }

for(i=0;i<600;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;i<600;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_B[MAXR][MAXC]={{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;i for(j=0;j    printf("\\n");

  }

  printf("按任一健继续...\\n");

  getch();

while(c for(i=r;i=0.0001)break;

if(i if(i!=r)for(j=0;j    x0=M_B[i][j];

    M_B[i][j]=M_B[r][j];

    M_B[r][j]=x0;

      }

      x0=M_B[r][c];

for(j=0;j for(i=r+1;i         x0=M_B[i][c];

for(j=0;j      }

      r++;

for(i=0;i for(j=0;j         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;i      x0=M_B[i][c];

for(j=0;j    }

for(i=0;i for(j=0;j      printf("\\n");

    }

    printf("按任一健继续...\\n");

    getch();

  }

  printf("完毕,按任一健退出...\\n");

  getch();

}

void f2(int n)

{

  int i,j;

for(i=0;i    if(i!=j)M_B[i][n+j]=0;

    else M_B[i][n+j]=1;

  }

  f1(n,2*n);

}

void main(void){  f1(3,5);  }

程序LAB6_3.C 模拟产生服从N(120,400)分布的随机变量800次, 并画出统计直方图

#include

#include

#include

#include

int Np[4000]={/*下列数据为标准正态分布函数值×10000, 由程序LAB2_3.C 产生*/

    5000,5003,5007,5011,5015,5019,5023,5027,5031,5035,…};

int f(int p)

{

  int i;

for(i=0;i<4000;i++)if(Np[i]>=p)break;

  return i;

}

void main(void)

{

  int i,j,x;

  float F,X[800];

  randomize();

  printf("\\n\\n");

for(i=0;i<800;i++){

    x=random(10000);

if(x<5000)F=-f(10000-x);

    else F=f(x);

    F*=0.001;

    X[i]=120+20*F;

    printf("%8.2f",X[i]);

  }

  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;F<150;F+=3){

    x=0;

for(i=0;i<800;i++)if(X[i]>=F&&X[i]    rectangle(j,400-x,j+30,400);

    j+=30;

  }

  getch();

  closegraph();

}

程序2000A01.C 统计文件2000A1.txt中ATCG的个数

#include

void main(void)

{

  int n=-1,ATCG[40][4]={0,0};

  char c;

  FILE *fp;

  fp=fopen("2000A1.txt

  while((c=fgetc(fp))!=EOF){

    if(c=='.')n;

    else if(c=='a')ATCG[n][0];

    else if(c=='t')ATCG[n][1];

    else if(c=='c')ATCG[n][2];

    else if(c=='g')ATCG[n][3];

  }

  fclose(fp);

for(n=0;n<40;n)printf("{%d,%d,%d,%d},\\n",ATCG[n][0],

ATCG[n][1],ATCG[n][2],ATCG[n][3]);

  getch();

}