3.2.1. Равномерное распределение
)
,
( b
a
U
( )
⎪⎩
⎪
⎨
⎧
≥
≤
<
∈
<
<
−
=
,
,
0
,
,
,
,
1
b
x
a
x
b
a
R
b
a
b
x
a
a
b
x
f
ξ
Алгоритм 1.
α
ξ
)
(
a
b
a
−
+
=
3.2.2. Нормальное (гауссовское) распределение
)
,
(
2
σ
a
N
0
,
,
2 1
)
(
2 2
2
)
(
2 2
>
∈
=
−
−
σ
πσ
σ
ξ
R
a
е
x
f
а
х
Алгоритм 1. Зарезервирована константа
π
2
=
c
. 1)
1
ln
2
α
−
=
r
; 2)
2
α
ϕ
c
=
;
3)
ϕ
ξ
cos
1
r
=
,
ϕ
ξ
sin
2
r
=
; 4)
1 1
σξ
ξ
+
= a
,
2 2
σξ
ξ
+
= a
Алгоритм 2. 1)
1 2
1 1
−
=
α
v
,
1 2
2 2
−
=
α
v
; 2)
2 2
2 1
v
v
s
+
=
; 3) Если
1
≥
s
, вер- нуться к п. 1; 4)
s
s
r
/
)
ln
2
(
−
=
; 5)
r
v
1 1
=
ξ
,
r
v
2 2
=
ξ
; 6)
1 1
σξ
ξ
+
= a
,
2 2
σξ
ξ
+
= a
В качестве случайного числа можно взять любое из чисел
2 1
,
ξ
ξ
3.2.3. Экспоненциальное распределение
)
(
λ
E
( )
⎪⎩
⎪
⎨
⎧
<
>
≥
=
−
−
−
0
,
0
,
0
,
0
,
1 1
x
x
e
x
f
x
λ
λ
λ
ξ
Алгоритм 1.
α
λ
ξ
ln
−
=
3.2.4. Распределение
2
χ
с
k степенями свободы
)
(
1
k
H
( )
⎪
⎪
⎩
⎪
⎪
⎨
⎧
≤
∈
>
⎟
⎠
⎞
⎜
⎝
⎛
⎟
⎠
⎞
⎜
⎝
⎛
Γ
=
−
−
0
,
0
,
,
0
,
2 2
2 1
2 1
2
x
Z
k
x
e
x
k
x
f
x
k
ξ
Алгоритм 1.
∑
=
=
k
i
i
u
1 2
ξ
, где
k
u
u
,...,
1
– независимые случайные величины с распределением
)
1
,
0
(
N
3.2.5. Распределение Стьюдента с
k степенями свободы
)
(
1
k
T
Z
k
k
x
k
k
k
x
f
k
∈
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
+
⎟
⎠
⎞
⎜
⎝
⎛
Γ
⎟
⎠
⎞
⎜
⎝
⎛ +
Γ
=
+
−
,
1 2
2 1
)
(
2 1
2
π
ξ
Алгоритм 1.
k
v
u
/
=
ξ
, где u и v – независимые случайные величины, при- чем
)
1
,
0
(
N
u
∈
,
)
(
1
k
H
v
∈
3.2.6. Распределение Фишера с
k
m, степенями свободы
)
,
(
1
k
m
F
( )
(
)
⎪
⎪
⎩
⎪
⎪
⎨
⎧
≤
∈
>
+
⎟
⎠
⎞
⎜
⎝
⎛
Γ
⎟
⎠
⎞
⎜
⎝
⎛
Γ
⎟
⎠
⎞
⎜
⎝
⎛ +
Γ
=
+
−
0
,
0
,
,
,
0
,
2 2
2 2
1 2
2 2
x
Z
m
k
x
mx
k
x
m
k
m
k
m
k
x
f
m
k
m
m
k
ξ
Алгоритм 1.
k
w
m
v
/
/
=
ξ
, где v и w – независимые случайные величины с рас- пределениями
)
(
1
m
H
и
)
(
1
k
H
соответственно.
3.2.7. Одномерное распределение Уишарта
)
,
(
2 1
σ
k
W
( )
⎪
⎪
⎩
⎪
⎪
⎨
⎧
≤
∈
>
>
⎟
⎠
⎞
⎜
⎝
⎛
Γ
=
−
−
0
,
0
,
,
0
,
0
,
2 2
1 2
2 1
2 2
2
x
Z
k
x
e
x
k
x
f
x
k
k
k
σ
σ
σ
ξ
Алгоритм 1.
∑
=
=
k
i
i
u
1 2
ξ
, где
k
u
u
,...,
1
– независимые случайные величины с распределением
)
,
0
(
2
σ
N
3.3. Средства Matlab для моделирования одномерных случайных чисел
3.3.1.Создание m-файлов-функций
В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пункт
Файл, затем пункты
Создать и
М-файл. После раскрытия окна редактора-отладчика необ- ходимо набрать нужные команды программы, отредактировать их и сохранить полученный файл под именем
f_name с помощью пунктов меню
Файл, Сохра-нить как… редактора-отладчика.
3.3.2. Управляющие структуры языка программирования системы Matlab 3.3.2.1. Диалоговый ввод-вывод
Disp(X) отображает массив, не печатая имя массива. Если X – строка, то ото- бражается текст.
Пример
x=[1 2 3]; disp(x)
1 2 3 disp('квадрат второго элемента=') квадрат второго элемента= disp(x(2)^2)
4
R=INPUT('Сколько яблок?')
дает пользователю приглашение в текстовой строке и затем ожидает ввода с клавиатуры. Может быть введено любое Matlab- выражение, которое вычисляется, используя переменные в текущей рабочей области, и результат возвращается в R. Если пользователь нажимает клавишу возврата каретки, ничего не вводя, то вводится пустая матрица.
R=INPUT('введите ваше имя’ , 's') дает приглашение в текстовой строке и ожидает ввода символьной строки. Напечатанный текст не вычисляется; сим- волы просто возвращаются как Matlab-строка.
Пример
r=input('введите угол в радианах') введите угол в радианах 2*pi r =
6.2832 r=input('введите ваше имя','s') введите ваше имя 2*pi r =
2*pi
3.3.2.2.
Циклы типа for-end
Циклы типа for-end обычно используются для организации вычислений с
заданным числом повторений цикла. Конструкция такого цикла имеет вид:
for var=выражение
Инструкция,…, Инструкция
end Выражение чаще всего записывается в виде b:s:e, где b – начальное значение переменной цикла var, s – приращение (шаг) этой переменной и e –
конечное значение управляющей переменной, при достижении которого цикл завершает- ся. Возможна запись выражения в виде b:e, в этом случае s=1. Список выпол- няемых в цикле инструкций завершается оператором
end.
Для досрочного выполнения цикла можно использовать оператор
break. Как только этот оператор встречается в программе, цикл прерывается.
Возможно использование цикла в цикле.
Пример 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
3.3.2.3. Циклы типа while_end
while Условие
Инструкции
end
Цикл типа while выполняется до тех пор, пока выполняется Условие. Для прекращения выполнения цикла можно использовать оператор break.
Пример
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
3.3.2.4.
Условный оператор 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
В качестве Оператора_отношения используются следующие логические операторы: ==, <, >, <=, >=, =. Двойные символы не имеют между собой про- белов.
Пример
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
3.3.2.5.
Переключатель 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.
Пример
Пусть существует 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
Ошибка в данных
3.3.2.6.
Создание паузы в вычислениях
Для остановки программы используется оператор pause в следующих фор- мах:
pause –
останавливает вычисления до нажатия любой клавиши.
pause(N) –
останавливает вычисления на N секунд.
pause on
– включает режим создания пауз.
pause off
– выключает режим создания пауз.
3.3.3. Стандартные функции Matlab для моделирования одномерных
случайных чисел
y=chi2rnd(k) –
2
χ
-распределение.
y=exprnd (lambda) – экспоненциальное распределение.
y=frnd (m,k) – распределение Фишера.
y=gamrnd (a,b) – гамма-распределение.
y=normrnd (a,sigma) – нормальное распределение.
y=trnd (k) – распределение Стьюдента.
y=unifrnd (a,b) – равномерное распределение.
y=rand(m,k)
– моделирует
)
(
k
m
×
-матрицу со случайными данными, вы- бранными из равномерного распределения в интервале (0,1).
r=unidrnd(k)
– возвращает матрицу случайных чисел, выбранных из набора
}
,...,
2
,
1
{
k . Размер r является размером k .
r=unidrnd(k,mm,nn) – возвращает
)
(
nn
mm
×
-матрицу случайных чисел, вы- бранных из набора
}
,...,
2
,
1
{
k .
3.4. Порядок выполнения работы
3.4.1.
Выполнить моделирование случайных чисел с указанными в п. 3.2 распределениями. Для каждого распределения вывести по 100 случайных чи- сел, используя собственную программу, реализующую предложенный алго- ритм, и стандартную программу Matlab. Собственные программы оформить в виде m-файлов-функций. Случайные числа вывести в виде точек на действи- тельной прямой.
3.4.2.
Для каждой выборки вычислить с помощью функции
function [xmean,s2,s3,s4,xmin,xmax,wtsum,wt,iwt,ifail]=g01aaf(x<,wt,iwt,ifail>)
выборочные среднее xmean, среднее квадратичное отклонение s2, коэффициент асимметрии s3, коэффициент эксцесса s4, минимальное значение выборки
xmin
, максимальное значение выборки xmax, сумму весов wtsum по данным
x
1
,x
2
,…,x
n
, помещенным в векторе x и имеющим соответствующие веса
w
1
,w
2
,…,w
n
, помещенные в векторе wt. Если присваивания весов не требуется, то параметр wt не указывается, при этом веса устанавливаются равными 1. Па- раметр iwt=0.
Пример n=4; iwt = 0; x = [193.0 215.0 112.0 161.0];
[xmean,s2,s3,s4,xmin,xmax,wtsum,wt,iwt,ifail] = g01aaf(x); xmean,s2,s3,s4,xmin,xmax,wtsum,ifail
ЛАБОРАТОРНАЯ РАБОТА № 4. МОДЕЛИРОВАНИЕ МНОГОМЕРНЫХ
СЛУЧАЙНЫХ ЧИСЕЛ
4.1. Цель работы
4.1.1.
Изучение методов моделирования многомерных случайных чисел.
4.1.2.
Приобретение навыков моделирования многомерных случайных чи- сел в системе Matlab.
4.2. Теоретические положения
4.2.1. Моделирование случайных чисел с многомерным нормальным
(гауссовским) распределением
Пусть требуется моделировать значения m -мерного случайного вектора
)
,...,
(
1
m
ξ
ξ
ξ
=
, распределенного по нормальному закону
)
,
( R
A
N
с математиче- ским ожиданием A и ковариационной матрицей R,
)
(
i
a
A
=
,
)
(
, j
i
R
σ
=
,
m
j
i
,
1
,
=
Обозначим
)
,...,
(
1
m
η
η
η
=
стандартный гауссовский m-мерный случайный вектор, то есть случайный вектор, распределенный по нормальному закону
)
,
0
( I
N
, где I – единичная матрица. Методы моделирования базируются на сле- дующем результате.
Теорема. Пусть
)
(
, j
i
c
C
=
– действительная
)
(
m
m
×
-матрица, являющаяся решением матричного уравнения
R
CC
T
= . (4.1)
Тогда случайный вектор
ξ
, являющийся линейным преобразованием
η
,
A
C
+
=
η
ξ
, (4.2) имеет нормальное распределение
)
,
( R
A
N
Моделирование
η
легко осуществляется. Действительно, компоненты этого вектора
m
η
η
,...,
1
некоррелированы, следовательно, и независимы, распределе- ние отдельной компоненты
i
η
– стандартное нормальное
)
1
,
0
(
N
. Поэтому мо- делирование
η
можно выполнить m -кратным обращением к функции модели- рования случайного числа с одномерным стандартным нормальным распреде- лением
)
1
,
0
(
N
Различные методы моделирования
ξ
, известные в литературе, отличаются лишь способом построения матрицы C.
Один из методов использует неединственность решения (4.1) и требует, чтобы C была нижней треугольной матрицей: 0
,
=
j
i
с
, если
i
j
>
. Ненулевые элементы
j
i
с
,
определяются рекуррентно. Действительно, в силу (4.2)
1 1
1
,
1 1
a
c
+
=
η
ξ
. Выражение (4.1) дает соотношение
1
,
1 2
1
,
1
σ
=
c
. Следовательно,
1
,
1 1
,
1
σ
=
c
.Теперь из (4.2) получим
2 2
2
,
2 1
1
,
2 2
a
c
c
+
+
=
η
η
ξ
, а из (2.1) имеем
,
2
,
1 1
,
1 1
,
2 2
,
2 2
2
,
2 2
1
,
2
σ
σ
=
=
+
c
c
c
c
Отсюда
1
,
1 2
,
1 1
,
2
c
c
σ
=
,
2 1
,
2 2
,
2 2
,
2
c
c
−
=
σ
, и т.д. Справедлива общая рекуррентная формула
∑
∑
−
=
−
=
−
−
=
1 1
2
,
,
1 1
,
,
,
,
/
)
(
j
j
j
j
j
j
i
j
i
j
i
c
c
c
c
ν
ν
ν
ν
ν
σ
σ
. (4.3)
Здесь
0
)
(
0 1
∆
=
=
⋅
∑
ν
и вычисления по рекуррентной формуле (4.3) осуществляются по строкам мат-
рицы
С , то есть в следующем порядке: ,
1
,
1
c,
,
2
,
2 1
,
2
cc,
,
,
3
,
3 2
,
3 1
,
3
cccmmcc,
1
,
4
,...,
Таким образом, моделирующий алгоритм определяется формулами (4.3),
(4.2).
Скачать 1.26 Mb.