Пример.
Вычислить интеграл
Методом Монте-Карло .
Решение.
Выберем случайный вектор
с плотностью распределения
,
равной
=
.
Для
выполняются условия
>0,
,
.
Ее можно рассматривать как совместную
плотность распределения двух независимых
экспоненциально распределенных
случайных величин
и
с плотностями распределения
и
Тогда
Пусть в результате моделирования
независимых компонент случайного
вектора
получена
случайная выборка
.
Согласно методу Монте-Карло в качестве
приближенного значения интеграла
следует использовать статистику
Текст программы:
#include <fstream.h>
#include <math.h>
#include <stdlib.h>
#include <limits.h>
#include <time.h>
// basic variate
double bv(double*a, double*b, double*c,int m,int n){
int i, tmp;
double result = c;
for =(i =0;i<n,i ++)
resalt +=b[i]*a[i];
tmp=(int)(resalt/m);
resalt -=((double)tmp)*((double)m);
for(i=1;I<n;i++)
a[i-1]=a[i];
a[n-1]=resalt;
return resalt/m;
}
// exponential variate p(x )=lexp(-lx),x>=0
double ev(double 1,double a) {
return –log(1-a)/1;
}
double avg(double*a,int n){
double s=0;
for(int i=0;i<n;i++)
s+=*(a+i);
return s/n;
}
double var(double*a,int n){
double s=0,s=2;
for(int i=0;i<n;i++){
s+=a[i];
s2+=a[i]*a[i];
}
return(s2-s*s/n)/(n-1);
}
void main() {
int i,j;
const int n=1000;
double a[2][n];//two pseudorandom sequences
double b[2][n];//coefficients
double c[2]; //increment
int m[2] //modules
//initial values
srand((unsigned) time (NULL) );
for(j=0;j<2;j++){
m[j] =RAND_MAX-rand();
for(i=0;i<n;i++){
a[j][i]=rand()%m[j];
b[j][i]=rand();
}
c[j]=rand();
}
//computation
double tmp;
double*A =new double[n];
double*x =new double[n];
double*y =new double[n];
ofstream out(“integral.txt”);
for(j=1;j<=5;j++){
for(i=0;i<n;i++){
do{
do{
x[i]=ev(2,dv(a[0],b[0],c[0],m[0],n));
y[i]=ev(0.5,dv(a[1],b[1],c[1],m[1],n));
}while(x[i]-y[i]==0);
tmp= 1+log(pow(y[i]-x[i],2));
}while(tmp<0);
A[i]=sqrt(tmp);
}
out<<avg(A,n)<<”\n”;
out<<var(A,n)<<”\n”;
out<<0.6745*sqrt(var(A,n)/n)<<”\n\n”;
}
}
Некоторые статистические результаты моделирования:
Запустим программу по 5 раз для различных значений n.
№ таблицы |
№ эксперимента
|
n |
Ответ |
Дисперсия |
Вероятная ошибка |
1 |
1 |
10 |
1.48915 |
0.268116 |
0.110444 |
|
2 |
|
1.37017 |
0.446164 |
0.142472 |
|
3 |
|
1.33855 |
0.097724 |
0.0656633 |
|
4 |
|
1.49458 |
0.13842 |
0.0793562 |
|
5 |
|
1.45071 |
0.334751 |
0.123408 |
2 |
6 |
20 |
1.4981 |
0.322398 |
0.0856373 |
|
7 |
|
1.21044 |
0.326385 |
0.0861652 |
|
8 |
|
1.55455 |
0.21716 |
0.070284 |
|
9 |
|
1.3085 |
0.366371 |
0.0912909 |
|
10 |
|
1.47881 |
0.150203 |
0.0584529 |
3 |
11 |
50 |
1.43942 |
0.268284 |
0.0494077 |
|
12 |
|
1.47837 |
0.217818 |
0.0445189 |
|
13 |
|
1.40335 |
0.258855 |
0.045317 |
|
14 |
|
1.43447 |
0.295186 |
0.0518256 |
|
15 |
|
1.42966 |
0.288469 |
0.0512326 |
4 |
16 |
100 |
1.39962 |
0.219702 |
0.0316154 |
|
17 |
|
1.43421 |
0.254309 |
0.0340144 |
|
18 |
|
1.42465 |
0.236301 |
0.032788 |
|
19 |
|
1.38404 |
0.253826 |
0.0339821 |
|
20 |
|
1.52162 |
0.257893 |
0.0342533 |
5 |
21 |
500 |
1.41239 |
0.272168 |
0.0157368 |
|
22 |
|
1.42945 |
0.271058 |
0.0157046 |
|
23 |
|
1.41602 |
0.259105 |
0.0153545 |
|
24 |
|
1.43384 |
0.243818 |
0.0148946 |
|
25 |
|
1.40812 |
0.246367 |
0.0149723 |
6 |
26 |
1000 |
1.42914 |
0.259182 |
0.00343388 |
|
27 |
|
1.42468 |
0.260844 |
0.00344486 |
|
28 |
|
1.42849 |
0.255745 |
0.00341103 |
|
29 |
|
1.42793 |
0.262582 |
0.00345633 |
|
30 |
|
1.42642 |
0.257236 |
0.00342295 |
Вычислить интеграл аналитически не удается. Полученное с помощью программы Mathematica4.1 точное значение интеграла равно 1.4221. Как видим, полученные в результате выполнения программы значения интеграла близки к точному, однако точность в данном случае невысока. Точность возрастает с увеличением объема выборки. Вероятная ошибка уменьшается.