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

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

Имитационное моделирование

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

www.domath.ru

Плотность и функция распределений

>    restart;

>    with(plots):#подключаем пакеты: plots пакет для построения графических объектов,
with(stats);#stats пакет функций по работе со стат. данными и
with(statplots);# пакет statplots пакет для построения стат. графики

>    N:=15:#зададим параметр

>   

# Зададим отрезок [a,b]

a:=-5+0.75/N:
b:=(N+19.25)-(N+1)^2/100:

>    #плотность и функцию распределения зададим как кусочно-непрер. фц.

>    f:=x->piecewise(a<x and x<b,2*(x-a)/(b-a)^2,a>x and x>b,0);
F:=x->piecewise(x<a,0,a<x and x<b,(x-a)^2/(b-a)^2,x>b,1);

Графики функций и плотности распределений

>    #зададим массив,где перечислены характеристики случайной величины
#необходимо для отображения полученных вычислений
mas:=Array(1..11,[

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

]):

Вычисления

>    #процедура, которая вычисляет момент k-го порядка непрерывной случайной величины (н.с.в.)

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


#зададим процедуру, зав. от двух параметров (a,n)
#данная процедура необходима для суммирования элементов массива a, заданной размерности n
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:

>    #вычисление математического ожидания

>   

#момент первого порядка есть мат. ожидание

M:=proc_un(f(x),0,1):
#вычисление среднеквадратического отклонения
sigma:=sqrt(proc_un(f(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(x),0,i));
                        print(mas[8] ,i, равен,proc_un(f(x),M,i));

end do;

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

>    #получение выборки объёма 900 равномерно распределённой случайной величины на интервале (0,1)

>    M:=900:
rany:=[stats[random,uniform[0,1]](M)]:

histogram(rany):

>    ran:=sort(rany):
# функция сортировки (от меньшего к большему)
#вычисление реализаций xx[i] из условия, что величина распределена равномерна на (a,b)

>    for i from 1 to M do
                              xx[i]:=a+(b-a)*sqrt(ran[i]);
end do:

   
    
             X_min:=xx[1]:#минимальное имаксимальное значения
                 X_max:=xx[M]:


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

>    t:=X_max-X_min:
                  l:=trunc((1+3.32*log[10](M))):#длина
                  h:=t/l:#шаг
X[1]:=X_min:
X[l+1]:=X_max:

>    #формирование нового масива с шагом

>    for i from 1 to l-1 do
                               X[i+1]:=X[i]+h;
end do:

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

>   

#Двойной цикл, который считывает число попавших в полуинтервал [X[j],X[j+1]]
for j from 1 to l do
                 m[j]:=0:# перед началом второго цикла эти параметры должны быть нулевыми

                 for i from 1 to M do

                               if (X[j]<=xx[i]) and (xx[i]<X[j+1])
                                                          then m[j]:=m[j]+1;
#внутри второго цикла происходит накапливание суммы m,
#если xx[i] попадает в полуинтервал
                               end if;
                    
                 end do:
#перед концом второго цикла подсчитаем абсолютные частоты, как отношение к величине M, что есть объём выборки
               row[j]:=m[j]/M;#частоты
               data[j]:=Weight(X[j]..X[j+1],row[j]);
 end do:

>    #формирование массива необх. для построения полигона

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

>    r:=plot(f(x),x=-5..30,thickness=5,color=red):

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

>    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,7), 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)
):#кусочно-непрерывная функция

>    plot([F(x),F1(x)],x=-15..X_max+3,
color=[green,red],
legend=["Theoretical","Empirical"],
style=[line,point]
);

>    proc_dis:=proc(k,m,t):
                             add((X[i]-m)^k*t[i],i=1..l);
end proc:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[Maple Plot]

[Maple Plot]

>    #вычисление основных выборочных характеристик

>    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], `равно`, 17.74966549

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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