Материал: Matlab. Практический подход. Самоучитель

Глава 8. Обработка данных

Рис. 8.24. Генерирование случайных чисел с показательным распределением

Использованный в вычислениях командный код приведен ниже (для удобства восприятия код снабжен комментариями):

%Количество тестов: n=10000;

%Параметр распределения: alpha=3;

%Равномерно распределенные случайные числа: xi=rand(1,n);

%Числа с показательным (экспоненциальным) распределением: eta=-log(1-xi)/alpha;

%Наибольшее из чисел:

Xmax=max(eta);

341

Самоучитель Matlab

%Количество рангов в гистограмме: k=100;

%Шаг дискретности:

h=Xmax/k;

% Аргумент для функции плотности распределения: x=0:h:Xmax;

%Создание гистограммы: hist(eta,x);

%Переход в режим удержания графики: hold on;

%Коэффициент масштабирования для плотности распределения: R=n*h;

%Значения масштабированной плотности распределения: p=Xmax*n/k*alpha*exp(-alpha*x);

%График функции (масштабированная плотность распределения): plot(x,p,'r-','LineWidth',2);

%Отображение сетки:

grid on;

% Отображение заголовка:

title('Показательное (экспоненциальное) распределение');

%Отображение легенды: legend('статистика','теория');

%Выход из режима удержания графики: hold off;

Переменная n определяет количество генерируемых случайных чисел. В принципе, чем это число больше, тем статистическое распределение ближе к теоретической кривой. Параметр распределения α записывается в переменную alpha командой alpha=3 (то есть в нашем случае рассматривается распределение с параметром α = 3 ). Командой xi=rand(1,n) создаем список из равномерно распределенных на интервале от 0 до 1 чисел. Числа с показательным (экспоненциальным) распределением создаем с помощью команды eta=-log(1-xi)/alpha.

Собственно, на этом процесс генерирования случайных чисел с нужным нам распределением завершен. Дальше мы пытаемся проверить, насколько сгенерированные числа (имеются в виду числа в переменной eta – именно эти числа мы будем использовать в дальнейшем) подчиняются экспоненциальному распределению. Делать это будем с помощью графиков. В частности, по набору случайных чисел построим гистограмму распределения случайных чисел по значениям и теоретическую кривую на основе функции плотности экспоненциального распределения.

Для построения графика функции плотности распределения нам понадобится наибольшее из вычисленных чисел (для верхней границы по

342

Глава 8. Обработка данных

аргументу при отображении функции). Наибольшее из случайных чисел вычисляем командой Xmax=max(eta). Переменная k определяет количество рангов в гистограмме (столбцов в диаграмме на один больше). Шаг дискретности по аргументу определяется в соответствии с количеством рангов гистограммы командой h=Xmax/k. Список аргументов для функции плотности распределения вычисляется как x=0:h:Xmax. Эти же значения будем использовать для определения рангов гистограммы. Командой hist(eta,x) создаем гистограмму. Здесь мы использовали встроенную функцию hist() для создания гистограмм. Первым аргументом передается список значений, для которых строится гистограмма, а второй аргумент определяет ранговые значения, для формирования столбцов гистограммы.

На заметку

Если строится гистограмма на основе значений {ηi } (первый аргумент функции hist()), i = 1,2,...,n с ранговой разбивкой в соответствии со списком

{xj } (второй аргумент функции hist()), j = 1,2,...,m , то в результате строится диаграмма, высота столбиков которой равняется количеству значений из списка {ηi }, которые попадают в диапазон между соседними точками xj и xj +1 (для всех возможных индексов). Функция hist() в качестве результата возвращает список со значениями высоты столбцов гистограммы.

Гистограмма дает нам представление о законе распределения случайных чисел. Нам этот закон необходимо сравнить с теоретической кривой для экспоненциального распределения. Соответствующая теоретическая зависимость строится на основе функции плотности распределения. Так, если случайные числа имеют плотность распределения p(x), а всего чисел n , то оценкой (точнее, математическим ожиданием) для количества чисел, попадающих в интервал значений от x до x +dx (при малых значениях dx ), будет значение np(x)dx . В нашем случае роль dx играет переменная h.

ния R=n*h. Значения "масштабированной" плотности распределения вычисляем командой p=Xmax*n/k*alpha*exp(-alpha*x). Для сравнения командой plot(x,p,'r-','LineWidth',2) строим график теоретической зависимости. Приведенный код сохраняем в файле с именем MyDist.m и в командном окне выполняем команду MyDist. Результаты графических построений представлены на рис. 8.25.

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

343

Самоучитель Matlab

Рис. 8.25. Распределение случайных чисел

Статистические вычисления

- Хотите доказательств?

- Нет, не хочу. Сейчас начнете нудеть, пальцы загибать: "Это раз, это два, это три...". А у меня времени нет. Мне нужно быстренько жирную точку на этом деле поставить.

К/ф "Статский советник"

Существует ряд задач, при решении которых (с разной степенью успешности) могут применяться статистические методы. В качестве примера рассмотрим следующую задачу. Есть набор статистических данных (случайные числа). Известно, что они распределены по нормальному закону с некоторым средним значением a и дисперсией σ2 . Эти параметры неизвестны, и их нужно определить (на основе статистических данных).

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

344

Глава 8. Обработка данных

На заметку

Моментом порядка k случайной величины ξ называется математическое ожидание от величины ξk (обозначается как Mξk ). Момент первого порядка называется математическим ожиданием случайной величины. Если случайная величина имеет плотность распределения p(x), то момент порядка k для этой

лить два параметра, то нас будут интересовать первые два момента, то есть

ские оценки необходимо сравнить с эмпирическими (вычисляемыми по статистическим данным).

На заметку

Эмпирическим моментом порядка k , вычисляемым на основе выборки {xi } n

объема n (то есть i = 1,2,...,n ), называется число mk = 1 ∑xik .

n i=1

Учитывая все приведенные выше обстоятельства, приходим к выводу, что

Документ с исполнительным кодом приведен на рис. 8.26.

Помимо непосредственно вычислений, этот код содержит также и команды для генерирования случайных чисел с параметрами (3, 4) (то есть a = 3 и σ2 = 4 ). Рассмотрим этот код:

n=100000;

x=3+2*randn(1,n);

m1=sum(x)/n;

a=m1;

m2=sum(x.^2)/n; sigma=sqrt(m2-a^2); disp('Параметры распределения:');

disp(['a=',num2str(a),' и sigma=',num2str(sigma)]);

345