Смесь гаусовских распределений

 

Плотность распределения НСВ ξ представляет собой смесь двух гаусовских распределений N ₁ (μ _1, σ ₁ ²) и N ₂ (μ ₂ σ ₂ ²) (модель засорений Тьюки-Хьюбера с уровнем засорений ε Осуществить n реализаций моделирования CB ξ и провести графический анализ результатов моделирования при следующих значениях параметров:

1) μ₁=μ₂=0;_, σ₂=ασ_1,α=1/2,2,4, σ₁=1; ε=0.05,0.45,n=20,50,100,500

2) μ₁=0,μ₂=6;_, σ₂=ασ_1,α=1/2,1,2, σ₁=1; ε=0.05,0.45,n=20,50,100,500

Глава 6. Метод Монте-Карло и его применения

Общая схема метода Монте-Карло

 

Метод Монте-Карло – это численный метод исследования математических моделей сложных систем, основанный на моделировании случайных элементов и последующем статистическом анализе результатов моделирования.

Метод Монте-Карло реализуется с помощью следующих шагов:

1. подбирают такую случайную величину ξ, что выполняются условия

                                   ,

где a-скалярная величина, приближенное значение, которое надо найти;

2. моделируют n независимых реализаций случайной величины ξ: ;

3. по случайной выборке , строят выборочную оценку

                                         

Точность вычислений характеризуют либо с помощью доверительного интервала, определяемого интервалом

                                  

вероятность выполнения которого приблизительно равна β; (например при ), либо с помощью так называемой вероятной ошибки метода Монте-Карло:

                                     

где численный множитель 0,6745 является решением уравнений

                     

т.е. одинаково вероятны ошибки, больше, чем  и ошибки, меньшие, чем

Вычисление определенного интеграла методом Монте-Карло

 

Рассмотрим задачу приближенного вычисления интеграла

где - подмножество из  При n=1 имеем определенный интеграл вида

Заметим, что схема вычислений как многомерных, так и одномерных интегралов, абсолютно аналогична.

Пусть η – произвольная случайная величина с плотностью распределения вероятностей

Предполагается только, что существуют моменты случайных величин, встречающиеся ниже. Рассмотрим случайную величину, являющуюся функциональным преобразованием случайной величины η.

 

Можно показать, что  Поэтому в качестве приближенного значения интеграла можно использовать статистическую оценку , построенную в выборке из n независимых случайных величин

 

Пример.

Вычислить интеграл  Методом Монте-Карло.

Решение.

Выберем случайный вектор  с плотностью распределения , равной = . Для  выполняются условия

>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. Как видим, полученные в результате выполнения программы значения интеграла близки к точному, однако точность в данном случае невысока. Точность возрастает с увеличением объема выборки. Вероятная ошибка уменьшается.