Моделирование случайных величин и статистическая обработка выборки в MATLAB

Ковач Ирина Петровна

(ГБОУ гимназия № 1528)

Алгоритмы моделирования случайных величин

Случайные числа с различными законами распределения обычно моде­ли­ру­ют­ся с по­мощью преобразований одного или нескольких независимых значений ба­зо­вой слу­чай­ной величины. Базовая случайная величина α – это случайная вели­чи­на с рас­пре­делением R(0;1) (равномерным распределением в интерва­ле (0;1)). В любой сис­те­ме программирования имеется стандартная программа моде­ли­ро­ва­ния базовой слу­чай­ной величины. Независимые слу­чай­ные величины с рас­пре­де­ле­ни­ем R(0,1) будем обоз­начать символами α₁, α₂, ... Рассмотрим алгоритмы мо­де­ли­ро­ва­ния случайных ве­ли­чин.

1. Равномерное распределение R(a, b):

.

Алгоритм моделирования:X = a + (b – a)α.

2. Нормальное (гауссово) распределение N(m, σ):

.

Алгоритм моделирования 1.

Зарезервирована константа c = 2π:

                1)      ;

                2)      ;

                3)      , ;

                4)      ,

Алгоритм моделирования 2.

                1)      , :

                2)      ;

                3)      Если s ≥ 1, вернуться к п. 1;

                4)      ;

                5)      , ;

                6)      , .

В качестве случайного числа можно взять любое из чисел X₁ и X₂.

3. Распределение хи-квадрат с k степенями свободы.

Алгоритм моделирования:

,

где  – независимые случайные величины.

4. Распределение Стьюдента с k степенями свободы.

Алгоритм моделирования:

,

где ,  – независимые случайные величины.

5. Распределение Фишера с k₁, k₂, степенями свободы.

Алгоритм моделирования:

,

где ,  - независимые случайные величины.

Средства MATLAB для моделирования случайных чисел

В MATLAB можно написать программу и сохранить ее в виде m-файла-сценария или m-файла-функции с целью последующего многократного выпол­нения. m-файл-функция является типичным объектом языка программирования системы MATLAB. Структура m-файла-функции с одним выходным пара­метром выглядит так.

function var=f_name(список параметров)

% Основной комментарий

%Дополнительный комментарий

Тело файла с любыми выражениями

var=выражение

Здесь переменная var – выходной параметр, f_name – имя функции.

Функция возвращает свое значение var и может использоваться в мате­ма­ти­че­ских выражениях в виде f_name(список параметров).

Все переменные в теле файла-функции являются локальными, то есть дей­ству­ют только в пределах тела функции, в отличие от файла-сценария, все пере­мен­ные которого являются глобальными.

Правила вывода комментариев те же, что и у файлов-сценариев.

Последняя конструкция (var=выражение) вводится, если требуется, чтобы функция возвращала результат вычислений. Если m-файл-функцию завершает стро­ка, в конце которой точка с запятой (;), то для возврата значения функции ис­поль­зуется программный оператор return.

Если выходных параметров больше одного, то структура модуля имеет сле­дую­щий вид.

function [var1,var2,…]=f_name(список параметров)

%Основной комментарий

%Дополнительный комментарий

Тело файла с любыми выражениями

var1=выражение

var2=выражение

……………………

Здесь var1, var2, … – имена переменных, которые являются выходными параметрами.

Такую функцию нельзя использовать в математических выражениях, по­сколь­ку она возвращает не один результат. Данная функция используется (вы­зы­ва­ет­ся) как отдельный элемент программы в виде:

[var1,var2,…]=f_name(список параметров).

Если такая функция используется в виде f_name(список пара­мет­ров), то возвращается значение только первого выходного параметра – переменной var1.

Если внутри функции целесообразно использовать глобальные перемен­ные, то их нужно объявить с помощью команды

global var1 var2 …

Начиная с версии 5.0 в функции системы MATLAB можно включать под­функ­ции. Они имеют такую же структуру, как и основная функция, и записы­ваются в теле основной функции.

Для создания и отладки m-файла-функции необходимо войти в редактор-от­лад­чик MATLAB, выбрав в меню командного окна MATLAB пункт Файл, затем пункты Соз­дать и m-файл. После раскрытия окна редактора-отладчика необходимо на­брать нужные команды программы, отредактировать их и сохранить полученный файл под именем f_name с помощью пунктов меню Файл, Сохранить как… ре­дак­то­ра-отладчика.

Управляющие структуры языка программирования системы MATLAB

Диалоговый ввод-вывод

Disp(X) отображает массив, не печатая имя массива. Если X – текстовая строка, то отображается текст.

Пример 1.

x=[1 2 3];

disp(x)

1 2 3

disp('квадрат второго элемента=')

квадрат второго элемента=

disp(x(2)^2)

4

R=INPUT('Сколько яблок?') дает пользователю приглашение в тек­сто­вой строке и затем ожидает ввода с клавиатуры. Может быть введено любое MATLAB-выражение, которое вычисляется с использованием переменных в текущей ра­бо­чей области, и результат которого возвращается в R. Если пользователь на­жи­ма­ет кла­ви­шу <Enter>, ничего не вводя, то вводится пустая матрица.

R=INPUT('введите ваше имя', 's') дает приглашение в текстовой стро­ке и ожидает ввода символьной строки. Напечатанный текст не вычисляется, а сим­во­лы просто возвращаются как MATLAB-строка.

Пример 2.

r=input('введите угол в радианах')

введите угол в радианах 2*pi

r =6.2832

r=input('введите ваше имя','s')

введите ваше имя 2*pi

r =2*pi

Циклы типа for-end

Циклы типа for-end обычно используются для организации вычислений с за­дан­ным числом повторений цикла. Конструкция такого цикла имеет следующий вид.

for var=выражение

Инструкция,…, Инструкция

end

Выражение чаще всего записывается в виде b:s:e, где b – начальное зна­че­ние переменной цикла var, s – приращение (шаг) этой переменной и e – конечное зна­чен­ие управляющей переменной, при достижении которого цикл за­вер­шается. Воз­мож­на запись выражения в виде b:e, в этом случае s = 1. Список вы­пол­няе­мых в цикле инструкций завершается оператором end. Для досрочного вы­пол­не­ния цик­ла можно использовать оператор break. Как только этот опе­ра­тор по­яв­ля­ет­ся в про­грам­ме, цикл прерывается. Возможно использование цикла в цик­ле.

Пример 3.

for i=1:3

for j=1:3

a(i,j)=i+j;

end

end

a

a =

2  3  4

3  4  5

4  5  6

Циклы типа while_end

while Условие

Инструкции

end

Цикл типа while выполняется до тех пор, пока выполняется Условие. Для прекращения выполнения цикла можно использовать оператор break.

Пример 4

x=1;i=1;

while x<=3

y(i)=x;

x=x+0.5;

i=i+1;

end

y

y =1.0000  1.5000  2.0000  2.5000  3.0000

Условный оператор if-elseif-else-end

Условный оператор if в общем виде записывается следующим образом.

if Условие

       Инструкции 1

elseif Условие

Инструкции 2

else

       Инструкции 3

end

Эта конструкция допускает несколько частных вариантов. Простейший из них – следующий.

if Условие

       Инструкции

end

Данный оператор работает следующим образом. Пока Условие возвра­ща­ет логическое значение 1 (то есть выполняется), выполняются Инструк­ции. Опера­тор end указывает на конец списка Инструкций. Инструкции в списке раз­де­ля­ют­ся запятыми или точками с запятыми. Если Условие возвращает логическое значение 0 (то есть не выполняется), то Инструкции также не выполняются. Еще один вариант.

if Условие

       Инструкции 1

else

       Инструкции 2

end

В этом варианте выполняются Инструкции 1, если выполняется Условие 1, или Инструкции 2 в противном случае.

Условие в операторе if записывается в виде:

Выражение_1 Оператор_отношения Выражение_2

В качестве Оператора_отношения используются следующие логи­че­ские операторы: ==, <, >, <=, >=, ~=. Двойные символы не имеют между собой пробелов.

Пример 5.

for i=1:3

       for j=1:3

             if i==j

             a(i,j)=2;

             elseif abs(i-j)==1

             a(i,j)=-1;

             else

             a(i,j)=0;

             end

       end

end

a

a =

2  -1  0

-1  2  -1

0  -1  2

Переключатель switch-case-otherwise-end

Для осуществления множественного выбора (или ветвления) использу­ет­ся конструкция с переключателем типа switch.

switch switch_Bыражение

case case_Bыражение

Список_ инструкций

case { case_Bыражение1, case_Bыражение2,… }

Список_ инструкций

…

otherwise,

Список_ инструкций

End

Выполняется первый оператор case, у которого case_Bыражение со­от­вет­ств­ует switch_Bыражению. Если ни одно из case_Bыражений не соот­вет­ствует switch_Bыражению, то выполняется список инструкций после опера­тора otherwise (если он существует). Выполняется только один case, а потом выполнение про­дол­жа­ет­ся с оператора после end.

Пример 6.

Пусть существует m-файл-сценарий swit.m.

switch month

       case{1,2,3}

       disp('Первый квартал')

       case{4,5,6}

       disp('Второй квартал')

       case{7,8,9}

       disp('Третий квартал')

       case{10,11,12}

       disp('Четвертый квартал')

       otherwise,

       disp('Ошибка в данных')

end

Эта программа в ответ на значения переменной month (номер месяца) опре­де­ля­ет но­мер квартала и выводит сообщение. Как это происходит, видно из следую­щей про­грам­мы.

month=3;

swit

Первый квартал

month=10;

swit

Четвертый квартал

month=13;

swit

Ошибка в данных

Создание паузы в вычислениях

Для остановки программы используется оператор pause в следующих формах:

pause – останавливает вычисления до нажатия любой клавиши.

pause(N) – останавливает вычисления на N секунд. pause on – включает режим создания пауз.

pause off – выключает режим создания пауз.

Стандартные функции MATLAB для моделирования одномерных случайных чисел

y=chi2rnd(k) – -распределение.

y=exprnd (lambda) – экспоненциальное распределение.

y=frnd (k1,k2) – распределение Фишера.

y=gamrnd (a,b) – гамма-распределение.

y=normrnd(m,sigma) – нормальное распределение.

y=trnd (k) – распределение Стьюдента.

y=unifrnd (a,b) – равномерное распределение.

y=rand(m,k) – моделирует ()-матрицу со случайными данными, вы­бран­ными из равномерного распределения в интервале (0,1).

r = unidrnd(k) – возвращает матрицу случайных чисел, выбранных из на­бора {1,2,..., k}. Размер r является размером k.

r = unidrnd(k,mm,nn) – возвращает ()-матрицу случайных чи­сел, выбранных из набора {1,2,..., k}.

Первичная обработка выборки

Эмпирическая функция распределения. Пусть имеется выборка x₁, ..., x_n из распределения F_X(x). Простейший взгляд на нее состоит в том, что числа  x₁, ..., x_n считаются возможными значениями некоторой дискретной случай­ной величины , причем вероятности этих значений одинаковы и равны 1/n. Ряд распределений этой случайной величины имеет следующий вид.

Таблица 1

Ряд распределений случайной величины

Эмпирической, или выборочной, функцией распределения  назы­вает­ся функция распределения дискретной случайной величины :

.

В соответствии с этим определением эмпирическая функция распределения задается формулой

,

где n(x) – количество выборочных значений, меньших x; n – объем выборки. Эмпирическую функцию распределения удобно строить с использованием вариа­ци­он­ного ряда. В этом случае она определяется формулой

.

В этой формуле x⁽ⁱ⁾ – i-й член вариационного ряда, i = 1, 2, ..., n–1. Эмпириче­ская функция распределения  представляет собой ступенча­тую функцию, по­сколь­ку это функция распределения дискретной случайной величины – см. ри­су­нок 1.

Рис. 1. Эмпирическая и теоретическая функции распределения

Эмпирическая функция распределения является состоятельной оценкой гене­раль­ной (теоретической) функции распределения , так как, согласно теореме Гливенко, имеет место сходимость по вероятности

 при .

Гистограмма. Гистограмма – это ступенчатая фигура, огибающая кото­рой сверху является графиком кусочно-постоянной функции , являю­щей­ся оцен­кой генеральной (теоретической) плотности распределения вероятностей  f_X(x). Ги­стограмма строится следующим образом. Промежуток , содер­жа­щий все выборочные значения , делится на некоторое количество l не­пе­ре­се­каю­щихся промежутков длины , j = 1, 2, ..., l и подсчитывается количество вы­бо­роч­ных значений n_j, попавших в j-й промежутков. Если на каждом промежутке как на основании построить прямоугольник высотой , где  – длина j-го промежутка, то мы получим фигуру, которая называется гистограммой. Пример гистограммы приведен на рисунке 2.

Рис. 2. Гистограмма и теоретическая плотность распределения

Гистограмма является состоятельной оценкой генеральной плотности рас­пре­де­ления вероятностей при увеличении объема выборки n и числа промежутков l, если при этом стремится к нулю максимальная из длин промежутков разбиения.

Есть два способа построения гистограммы.

1. Равноинтервальный способ. Выбирают количество интервалов l, а дли­ну Δ каждого интервала определяют по формуле .

2. Равновероятный способ. Выбирают количество выборочных значений m, по­пав­ших в каждый промежуток. Объем выборки должен быть кратен m.
Тогда число про­ме­жутков l =n/m и интервалы будут следующими: , , …,  где x(i) – i-й член вариационного ряда. При этом способе промежутки имеют раз­лич­ную длину – и границы проме­жут­ков попадают на выборочные зна­че­ния. При­ня­то считать, что граничное значение делится поровну между двумя ин­тер­валами, то есть 1/2 значения попадает в левый промежуток и 1/2 – в правый. Понятно, что при этом в крайний левый промежуток попадает m – 1/2  значений, в крайний правый m + 1/2 значений, а в средние промежутки – по m значений.

Опишем программы для получения оценок законов распределения и отоб­ра­же­ния их на экране.

Сортировка в MATLAB

y = sort(x) сортирует элементы вектора x в возрастающем порядке. Здесь x – исходный вектор, y – отсортированный вектор. Если x – матрица, функция y = sort(x) сортирует каждый столбец x в возрастающем порядке.

y = sort(x,d) сортирует матрицу x вдоль измерения d.

[y,i] = sort(x,d) возвращает также индексную матрицу i. Если x – вектор, то элементы индексной матрицы указывают номера элементов вектора y в ис­ходном векторе x (см. toolbox\matlab\datafun).

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

Лестничные графики в MATLAB

stairs(y) строит лестничный (ступенчатый) график по значениям эле­мен­тов вектора y.

stairs(x,y) строит лестничный график по значениям элементов векто­ра y в точках скачков, определенных в x. Значения x должны распо­ла­гаться в воз­рас­таю­щем порядке (см. toolbox\matlab\specgraph).

Функция stairs(x,y) используется для получения и графического отоб­ра­же­ния эмпирической функции распределения. В этом случае x – ва­риа­ционный ряд, а y(i) = i/n, i = 1, 2, ..., n.

Гистограммы в MATLAB

n=hist(y) распределяет элементы вектора y в 10 интервалов одинаковой длины:

и возвращает количество элементов, попавших в каждый интервал, в виде вектора n. Если y – матрица, то hist работает со столбцами.

n=hist(y,l), где l – скаляр, использует l интервалов одинаковой длины:

.

n=hist(y,x), где x – вектор, возвращает количество элементов вектора y, попавших в интервалы с центрами, заданными вектором x. Число интервалов в этом случае равно числу элементов вектора x.

[n,x]=hist(…) возвращает числа попаданий в интервалы (в векто­ре n), а также положения центров интервалов (в векторе x).

hist(…) строит гистограмму без возвращения параметров, то есть строит прямоугольники высотой h_j = n_j, где n_j  – число элементов, попавших в j-й интервал, j = 1, 2, ..., l  (см. toolbox\matlab\datafun).

Функция hist используется для получения и отображения гистограм­мы.

Этапы статистической обработки выборки в MATLAB

                 1.       Создать m-файлы-функции, реализующие алгоритмы моделирования слу­чай­ных величин, приведенные в пункте 2.2.

                 2.       Выполнить моделирование случайных чисел с различными законами рас­пре­делений. Для каждого распределения вывести по 100 случай­ных чисел , используя:

а) собственные m-файлы-функции;

б) стандартную программу MATLAB.

                 3.       Для каждого распределения вывести на экран в одно графическое ок­но гис­тограмму и генеральную плотность распределения вероятно­стей, а в другое графическое окно – эмпирическую функцию распре­де­ления и гене­ральную функцию распределения. Для вывода гене­раль­ных плот­но­стей рас­пре­деления вероятности и функций распре­де­ле­ния исполь­зо­вать про­грам­мы MATLAB. Для согласования масш­та­бов гистограммы и гене­ральной плотности распределения вероят­ностей необходимо генеральную плотность распределения умножить на коэффициент

.

                 4.       Исследовать сходимость эмпирических распределений к генераль­ным при увеличении объема выборки n.

                 5.       Для каждой выборки вычислить с помощью функции:

Function

[xmean,s2,s3,s4,xmin,xmax,wtsum,wt,iwt,ifail]=g01aaf(x<,wt,iwt,ifail>);

·         выборочные среднее xmean;

·         среднее квадратичное отклонение s2;

·         коэффициент асимметрии s3;

·         коэффициент эксцесса s4;

·         минимальное значение выборки xmin;

·         максимальное значение выборки xmax;

·         сумму весов wtsum по данным , помещенным в векторе x и имею­щим соответствующие веса w1,w2, … , wn, помещенные в векторе wt.

Если присваивания весов не требуется, то параметр wt не указывается, при этом веса устанавливаются равными 1. Параметр iwt = 0.

                 6.       Сравнить найденные выборочные характеристики с соответствующими па­ра­метрами моделируемых распределений при различных объемах вы­бо­рок n

Литература

                 1.       Аллавердиев А.М., Ревякин А.М. Теория вероятностей и математическая ста­ти­сти­ка: Методическое пособие. М.: МГИДА, 2004.

                 2.       Афифи А., Эйзен С. Статистический анализ. Подход с исполь­зо­ва­ни­ем ЭВМ. М.: Мир, 1982.

                 3.       Вуколов Э.А. Основы статистического анализа: Практикум по стати­сти­ческим методам и исследованию операций с использованием па­ке­тов STATISTICA и EXCEL: Учебное пособие. М.: ИНФРА-М., 2004.

                 4.       Ревякин А.М., Терещенко А.М., Третьяков В.А. Методическое посо­бие по математической статистике для индивидуальной работы студен­тов вечернего факультета. М.: МИЭТ, 1989.

                 5.       Тюрин Ю.Н. и др. Теория вероятностей и статистика. М.: Изд-во МЦНМО; АО Московский учебник, 2004.

                 6.       Тюрин Ю.Н. и др. Теория вероятностей и статистика: Методическое пособие для учителя. М.: Изд-во МЦНМО МНОО, 2005.

                 7.        Теория и практика статистических исследований / Под. ред. А.М. Ревякина и В.В. Костылева. - М.: МГАДА, 2007.- 354 с.