Дисциплина:
Раздел:
Выполнил:

Математическая статистика

Имитационное моделирование. Непрерывная случайная величина.

ЧГУ, матфак, 5 курс Осипов.П.Г

www.domath.ru

>    restart:
with(plots):#пакет графики
with(stats):#пакет для работы со статистическими даными
with(statplots):# пакет для построения различных графиков по стат. данным

>    #массив характеристик случайной величины,
#которые требуется вычислить

mas:=Array(1..11,[

математическое_ожидание_M[X],
дисперсия,
среднеквадратическое_отклонение,
коэффициент_вариаций,коэффициент_симметрий,
медиана,начальный_момент_порядка,
центральный_момент_порядка,коэффициент_ассиметрии,
Коэффициент_эксцесса,
дифференциальная_энтропия
]

):

>    #процедура суммирования
summary:=proc(a,n) local s::float, i::int;

s:=0:          
for i from 1 to n
                 do
                  s:=s+a[i];
end do;

end proc:

 

>    #плотность и функция нормального распределения f_ksi:=x->1/(sigma*sqrt(2*Pi))*exp(-(x-M)^2/(2*sigma^2));
F:=x->Int(1/(sigma*sqrt(2*Pi))*exp(-(t-M)^2/(2*sigma^2)),t=-infinity..x);

#процедура предназ. для вычисления характеристик данного распределения
proc_un:=proc(func,m,k);
                         int((x-m)^k*func,x=-infinity..+infinity);
end proc:

f_ksi := x -> 1/sigma/sqrt(2*Pi)*exp(-1/2*(x-M)^2/sigma^2)

F := x -> Int(1/sigma/sqrt(2*Pi)*exp(-1/2*(t-M)^2/sigma^2),t = -infinity .. x)

>    (sigma,M):=(4.15,26.5):#начальные данные основных характеристик

>    plot(f_ksi(x),x=10..40);#график плотности
teo:=plot(F(x),x=10..40,color=green):#график функции распределения

>    M:=proc_un(f_ksi(x),0,1):#начальный момент первого порядка есть мат. ожидание
sigma:=sqrt(proc_un(f_ksi(x),M,2)):#центральный момент второго порядка есть среднк. отклонение
#вычисление остальных характеристик и вывод на экран
print(mas[1],равно,M);
print(mas[2],равно,sigma^2);
print(mas[3],равно,sigma);
print(mas[4],равно,(sigma/M)*100,процентов);

#print(mas[5],равно,M/sqrt(proc_dis(3,proc_dis(1,0,p),p)));
for i from 1 to 6
                     do
                        print(mas[7] ,i, равен,proc_un(f_ksi(x),0,i));
                        print(mas[8] ,i, равен,proc_un(f_ksi(x),M,i));

end do;

print(mas[9],равен,proc_un(f_ksi(x),0,3)/sigma^3);
print(mas[10],равен,proc_un(f_ksi(x),0,4)/sigma^4-3);
#print(mas[11],равна,-evalf(int(f_ksi(x)*ln(f_ksi(x)),x=-infinity..+infinity)));

`математическое_ожидание_M`[X], `равно`, 26.50000000

`дисперсия`, `равно`, 17.22250000

`среднеквадратическое_отклонение`, `равно`, 4.150000000

`коэффициент_вариаций`, `равно`, 15.66037736, `процентов`

`начальный_момент_порядка`, 1, `равен`, 26.50000000

`центральный_момент_порядка`, 1, `равен`, 0.

`начальный_момент_порядка`, 2, `равен`, 719.4724999

`центральный_момент_порядка`, 2, `равен`, 17.22250000

`начальный_момент_порядка`, 3, `равен`, 19978.81375

`центральный_момент_порядка`, 3, `равен`, 0.

`начальный_момент_порядка`, 4, `равен`, 566611.9097

`центральный_момент_порядка`, 4, `равен`, 889.8435187

`начальный_момент_порядка`, 5, `равен`, 16391556.09

`центральный_момент_порядка`, 5, `равен`, 0.

`начальный_момент_порядка`, 6, `равен`, 483168604.4

`центральный_момент_порядка`, 6, `равен`, 76626.65000

`коэффициент_ассиметрии`, `равен`, 279.5280585

`Коэффициент_эксцесса`, `равен`, 1907.263651

>    #объём выборки

>    n:=900:
#выборка заданного размера равномерно/распредел случ. вел.
rany:=seq([stats[random,uniform[0,1]](12)],i=1..n):

>    #обработка данных согласно методу имитац. моделирования для норм. случ. величины

>    for i from 1 to n do
                        xx[i]:=sigma*(add(rany[i,j],j=1..12)-6)+M;

end do:

>    #сортировка данных

>    ran:=sort([seq(xx[i],i=1..n)]):

   
    
                X_min:=ran[1]:
                  X_max:=ran[n]:


>    #размах выборки

>    t:=X_max-X_min:
                  l:=trunc((1+3.32*log[10](n)));# длина выборки
                  h:=t/l:#шаг выборки
X[l+1]:=X_max:

>    for i from 1 to l do
                    X[i]:=X_min+(i-1)*h:
                    w[i]:=(X[i]+X[i+1])/2:
                
end do:

>    #ниже заложен механизм им. моделирования.

>    for i from 1 to l do
                       m[i]:=0;
end do:  

>    for j from 1 to n do
                
                               if (ran[j]>=X[1]) and (ran[j]<X[2])
                                                          then m[1]:=m[1]+1;
                               end if;

                 for i from 2 to l do

                               if (ran[j]>X[i]) and (ran[j]<X[i+1])
                                                          then m[i]:=m[i]+1;
                               end if;
                
                 row[i]:=m[i]/n;
                 data[i]:=Weight(X[i]..X[i+1],row[i]);                end do:
end do:

>    row[1]:=m[1]/n:#частоты
data[1]:=Weight(X[1]..X[2],row[1]):

>    data:=[
data[1],data[2],data[3],data[4],data[5],data[6],data[7],data[8],data[9]]:

>    r:=plot(f_ksi(x),x=10..40,color=red):

>    t:=histogram(data):#полигон
display(t,r);#вывод на экран теоретической плотности распределения и полигон частот

[Maple Plot]

>    #эмпирическая функция распределения
F1:=x->piecewise(
x<=X[1],0,
X[1]<x  and x<=X[2], summary(row,1), X[2]<x  and x<=X[3], summary(row,2),
X[3]<x  and x<=X[4], summary(row,3), X[4]<x  and x<=X[5], summary(row,4),
X[5]<x  and x<=X[6], summary(row,5), X[6]<x  and x<=X[7], summary(row,6),
X[7]<x  and x<=X[8], summary(row,8), X[8]<x  and x<=X[9], summary(row,8),
X[9]<x  and x<=X[10],summary(row,9), X[10]<x,summary(row,10)
):

>    emp:=plot(F1(x),x=10..40,color=red,style=point):#график эмпирической функции распределения
display(teo,emp);#вывод на экран графиков теор. и эмпири. функций

[Maple Plot]

>    #процедура для вычисления характеристик эмпирич. распределения
proc_dis:=proc(k,m,t);
                             add((w[i]-m)^k*t[i],i=1..l);
end proc:

>    #математическое ожидание

>    M:=proc_dis(1,0,row):
sigma:=sqrt( proc_dis(2,M,row) ):#среднеквадратическое отклонение

>    #вычисление остальных характеристик и вывод на экран
print(выборочное,mas[1],равно,M);
print(выборочная,mas[2],равна,sigma^2);
print(выборочное,mas[3],равно,sigma);
print(выборочный,mas[4],равен,(sigma^2/M)*100,процентов);

#print(mas[5],равно,M/sqrt(proc_dis(3,proc_dis(1,0,p),p)));
for i from 1 to 6
                     do
                        print(выборочный,mas[7] ,i, равен,proc_dis(i,0,row));
                        print(выборочный,mas[8] ,i, равен,proc_dis(i,M,row));

end do;

print(выборочный,mas[9],равен,proc_dis(3,0,row)/sigma^3);
print(выборочный,mas[10],равен,proc_dis(4,0,row)/sigma^4-3);
print(выборочная,mas[11],равна,-evalf(add(ln(row[i])*row[i],i=1..l)));

`выборочное`, `математическое_ожидание_M`[X], `равно`, 26.73107767

`выборочная`, `дисперсия`, `равна`, 17.89959692

`выборочное`, `среднеквадратическое_отклонение`, `равно`, 4.230791524

`выборочный`, `коэффициент_вариаций`, `равен`, 66.96174820, `процентов`

`выборочный`, `начальный_момент_порядка`, 1, `равен`, 26.73107767

`выборочный`, `центральный_момент_порядка`, 1, `равен`, .2970119252e-1

`выборочный`, `начальный_момент_порядка`, 2, `равен`, 733.2440551

`выборочный`, `центральный_момент_порядка`, 2, `равен`, 17.89959692

`выборочный`, `начальный_момент_порядка`, 3, `равен`, 20575.85725

`выборочный`, `центральный_момент_порядка`, 3, `равен`, -2.720584086

`выборочный`, `начальный_момент_порядка`, 4, `равен`, 589636.3801

`выборочный`, `центральный_момент_порядка`, 4, `равен`, 901.9022092

`выборочный`, `начальный_момент_порядка`, 5, `равен`, 17229123.61

`выборочный`, `центральный_момент_порядка`, 5, `равен`, -8.1777318

`выборочный`, `начальный_момент_порядка`, 6, `равен`, 512647362.3

`выборочный`, `центральный_момент_порядка`, 6, `равен`, 69044.41472

`выборочный`, `коэффициент_ассиметрии`, `равен`, 271.7021415

`выборочный`, `Коэффициент_эксцесса`, `равен`, 1837.338742

`выборочная`, `дифференциальная_энтропия`, `равна`, 1.903456930