Численные операции над случайными величинами и их приложения Текст научной статьи по специальности «Математика»

УДК 519.24

Численные операции над случайными величинами и их приложения

Борис С.Добронец*

Институт космических и информационных технологий, Сибирский федеральный университет, Киренского, 26, Красноярск, 660074,

Россия

Ольга А. Попова^

Сибирская автомобильно-дорожная академия, Мира, 5, Омск, 644080, Россия

Получена 12.11.2010, окончательный вариант 01.12.2010, принята к печати 10.01.2011 В 'работе рассмотрено применение численных операций над случайными величинами. Приведены описание и свойства основных операций над плотностями случайных величин. Представлены алгоритмы и примеры решений систем линейных алгебраических и нелинейных уравнений, примеры использования развиваемых методов в практике принятия экономических решений. Показано, что данный подход позволяет существенно повысить качество принимаемых решений и сократить объем вычислений.

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

Введение

При решении большого числа практических задач возникает необходимость нахождения не только самих решений, но и различных оценок, таких как оценки погрешности, границы множеств решений и т. п.

Существующая неопределенность информации отражается в данных. Можно выделить три типа "неопределенных" данных: случайные, нечеткие и интервальные. Случайные числа задаются некоторыми вероятностными распределениями их возможных значений, нечеткие данные задаются лингвистически сформулированными распределениями их возможных значений, интервальные данные задаются интервалами их возможных значений без указания какого-либо распределения возможных значений числа внутри заданного интервала [1].

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

* [email protected] [email protected] © Siberian Federal University. All rights reserved

обобщения интервальной математики направлены на определение более детальных характеристик получаемых множеств решений [5].

Важное направление — вероятностное представление входных данных, например, в виде гистограммных чисел и разработка численных операций над ними [6, 7]. Одна из первых работ в этой области [8].

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

Сравнение подходов интервальной математики и численных операций над плотностями вероятности случайных величин показывает, что практически значимые плотности вероятности занимают лишь небольшую часть полученных интервалов. Поэтому в таких случаях границы решений мало информативны. Часто при работе со случайными величинами ограничиваются вычислением лишь небольшого числа характеристик: математическое ожидание и дисперсия.

Одним из подходов учета случайного характера входных данных является метод Монте-Карло [9,10]. При всех его положительных качествах этот метод обладает рядом недостатков. Один из самых существенных — низкая скорость сходимости.

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

1. Гистограммные переменные

Обозначим через К множество {ж} случайных величин, заданных своими плотностями ве-

п» «

роятности рх, соответственно, К — пространство случайных векторов.

Наряду с общими представлениями случайных величин своими плотностями в виде непрерывных функций, будем рассматривать случайные величины, плотность распределения которых представляет гистограмму. Гистограмма Р — кусочно-постоянная функция — определяется сеткой = 0, ...,п}, на отрезке [хг_1 ,хг] гистограмма принимает постоянное значение рг.

Рассмотрим вопрос построения гистограммы Р по некоторой рх. Тогда значение рг на отрезке [хг-1 ,хг] определится как среднее:

2. Законы распределения функций случайных аргументов

Имеется непрерывная случайная величина х с плотностью распределения рх. Другая случайная величина у связана с нею функциональной зависимостью:

Плотность распределения величины у в случае монотонной функции I согласно [11] определяется следующим образом:

У = I (х).

ру (у) = I (I-1(у))1(1 -1)'(у)1,

где /_1 — функция, обратная к /. Однако пользоваться таким представлением в ряде случаев довольно затруднительно.

Рассмотрим построение гистограммы ру по известной гистограмме рх. Пусть рх задается сеткой {х®|г = 0,..., п}. Вычислим сетку для гистограммы ру {у®|г = 0,..., п} как у® = /(х®). Тогда вероятность попадания случайной величины х в отрезок [х®_1, х®] равна вероятности попадания случайной величины у в отрезок [у^_1, у^]. Следовательно, гистограмма ру будет определяться сеткой {у®|г = 0,..., п} и соответствующими значениями ру® = рх.

Рассмотрим общий случай, когда / не является монотонной (рис. 1). Пусть {у®|г = 0,..., п} — сетка для гистограммы ру. Фиксируем отрезок [у®_1, у®], найдем ру®. Пусть [а^-, ], 2 = 1,..., т такие, что [у®_1, у®] = /([а^-, Ьj]). Тогда вероятность ру® попадания случайной величины у в отрезок [у®_1,у®] определяется, как

т '

ру® = у рх(^к.

^ а

Рис. 1. Построение гистограммы функции случайного аргумента

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

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

^ = /(х, у).

Тогда закон распределения Р2 величины г [11]:

Р* (г) = Р (/(х, у) < г).

Пусть, как и выше, гистограммар2 определяется сеткой {¿¿|г = 0, ...,п}. Определим область П® = {(х,у)|^ < /(ж, у) < ¿¿+1}. Тогда имеет вид

р« = JJр(х,у)йхйу/(г®+1 - ¿¿). - 231 -

3. Численные операции

Пусть имеется система двух непрерывных случайных величин (ж1, Ж2) с плотностью распределения р(ж1,ж2). Известны аналитические формулы для определения плотности вероятности результатов арифметических действий над случайными величинами. Например, для нахождения плотности вероятности рХ1+Х2 суммы двух случайных величин Ж1 + Ж2 используется соотношение [12]

рХ1 +Х2 (ж) = J р(ж — = J р(«,ж — (1)

— то — то

для нахождения плотности вероятности рХ1/Х2 частного двух случайных величин Ж1/Ж2

то 0

РХ1/Х2 = J «р(ж-у,-у)й-у — J «р(«, (2)

0 —то

Плотность вероятности рХ1Х2 произведения двух случайных величин Ж1Ж2 [11] определяется соотношением

то 0

рХ1Х2 = У(1/«)р(ж/«, — J (1/«)р(«,ж/у)^. (3)

0 —то

Однако эти формулы не всегда удобны для численных расчетов.

Основные принципы разработки гистограммных операций продемонстрируем на примере операции сложения. Пусть г = жх + Ж2, и носители жх — [01,02], Ж2 — [61,62], р(жх,ж2) — плотность распределения вероятностей случайного вектора (Ж1, Ж2). Заметим, что прямоугольник [01,02] х [61,62] — носитель плотности распределения вероятностей р(ж1,ж2) (рис. 2) и плотность вероятности г отлична от нуля на интервале [01 + 61, 02 + 62]. Обозначим ¿¿, г = 0,1,..., п — точки деления этого интервала на п отрезков. Тогда вероятность попадания величины г в интервал [¿¿, £¿+1] определяется по формуле

п»

где П = {(ж1, ж2)|^ ^ ж1 + ж2 ^ 1} [8]. И окончательно р^ имеет вид

р« = р(ж1 ,ж2)йж1^ж2/(г4+1 — ¿¿).

Рассмотренный выше подход обобщается на случай большего числа переменных. Пусть требуется найти гистограмму р2 суммы

г = 01Ж1 + 02Ж2 + ... + о„Ж„,

и пусть р(ж1, ж2,..., жп) — плотность распределения вероятностей случайного вектора (Ж1, Ж2,..., жп). Тогда вероятность попадания г в интервал (¿¿,£¿+1) соответственно равна [12]

Р < г < £¿+1) = J ."У р(ж1, ж2,..., жп)йж1^ж2....йжп,

то

п

п

Рис. 2. Построение гистограммы суммы двух случайных величин

где П® = {(хьх2, ...,х„)< а^1 + а2х2 + ... + а„х„ < ¿¿+1}, имеет вид

р« = J р(х1,х2, ...,хп)йх1Йх2....йх„/(г®+1 - ¿¿). п»

Рассмотрим в качестве примера операцию тах(х, у). Вероятность того, что тах(х, у) < г определяется через функцию распределения Р:

р (¿) = | Рх(о#у Ру

Используя функцию распределения Р, можно построить гистограмму для тах(х, у) на сетке {¿¿}. Тогда р® = Р(¿¿+1) - Р(¿¿).

4. Решение уравнений

Рассмотрим системы линейных или нелинейных уравнений

/®(х, к) = 0, г = 1, .. ., п,

где х € Д" — случайный вектор решения, к € Дт — случайный вектор коэффициентов, р&(£ь £2,..., £т) — функция плотности вероятности вектора коэффициентов и К — носитель Множество решений X имеет вид

X = {х|/(х, к) = 0,к € К}.

Каждому х € X можно сопоставить подмножество коэффициентов Кх С K:

К = {к|/(х,к)=0.}

Пусть необходимо найти вероятность Р (Хо) попадания решения х в некоторое подмножество Хо С X. Сопоставим Хо множество коэффициентов Ко = {Кх|х € Хо}.

Тогда вероятность

Ко

Нахождение подмножеств КХ с К в общем случае не простая задача. Но для ряда частных случаев вполне реализуемая.

Рассмотрим задачу нахождения корня одномерного уравнения /(ж, к) = 0. Предположим, что корень локализован на отрезке [0,6]. Заметим, что ж будет представлять собой случайную величину и необходимо найти плотность распределения Ж этой случайной величины. Далее, мы будем предполагать, что можно с достаточной точностью для любого £ € [0,6] вычислять плотность распределения ф2 случайной величины f (г, к).

Тогда Р(г) есть вероятность, что корень лежит левее (правее) точки г:

о

Р(*)=/ Ф*

—то

Действительно, перепишем уравнение /(ж, к) = 0 в виде ф(к) — ж = 0. Тогда ж = ф(к) и ф* есть плотность распределения ф(к), или, что то же самое, плотность распределения ж, сдвинутая на г.

Таким образом, плотность распределения корня ж = Р'(г). В тех случаях, когда Р(г) — не аналитическая функция и явное вычисление производной затруднено, можно вычислять производную численным способом. Для построения гистограммы д с узлами {¿¿}, приближающей плотность распределения ж, достаточно вычислить Р(¿¿). Тогда значение гистограммы g¿ на отрезке [¿¿_1,2^] определяется соотношением

Р (¿¿) — Р (¿¿_1) g¿ =---.

г ¿¿_1

В качестве примера рассмотрим решение простейшего квадратного уравнения аж2 — Ь = 0, где а, Ь — случайные величины с равномерным законом распределения, заданные, соответственно, на отрезках [1, 2], [2,4]. Несложно убедиться, что ж задана на отрезке [1, 2]. Таким образом, для Уж € [1, 2]: /(ж, а, Ь) = аж2 — Ь.

Заметим, аж2 — равномерная случайная величина, заданная на отрезке [ж2, 2ж2], тогда для вычисления / необходимо найти разность двух независимых равномерных случайных величин. Для построения гистограммы Ж вычислим функцию распределения РХ в точках сетки {ж¿ = 1 + ¿Н, Н = 1/п, г = 0,1, ...,п}. Далее по ней построим гистограмму для корня квадратного уравнения. На рис. 3 приведена гистограмма ж корня квадратного уравнения при п = 100.

Рассмотрим решение систем линейных алгебраических уравнений

Аж = 6.

Для простоты изложения предположим, что случайный вектор Ь состоит из независимых компонент Ь1, Ь2, каждая из которых случайная величина, равномерно распределенная на отрезке [0,1]. Тогда носитель плотности вероятности вектора Ь — квадрат [0,1]2. Пусть матрица А имеет вид

А=( —2

Построим функцию распределения случайного вектора Ж. Заметим вероятность того, что ж1 < Г1 и ж2 < Г2 равна площади, отсекаемой от квадрата [0,1]2 прямыми, проходящими через точку 6о = Аг, где г = (г1, Г2)т с направляющими векторами /1 = (2, —1) /2 = ( — 1, 2).

2 "

1 "

Рис. 4. Вычисление плотности вероятности

На рис. 4 в левой части серым цветом закрашено X — множество решений системы. Для вычисления вероятности попадания в квадрат к отобразим его с помощью матрицы А на вектор правой части Ь (правая часть рисунка). Квадрат к отображается в ромб Д. Таким образом, вероятность попадания решения в квадрат к равняется интегралу от плотности вероятности вектора Ь по области пересечения ромба Д и вектора правой части Ь.

Непосредственными вычислениями нетрудно получить плотность вероятности случайного вектора х. В левой части рис. 5 приведено кусочно-постоянное с шагом 0,1 приближение совместной плотности вероятности вектора х. Сплошной линией приведена граница множества решений исходной системы. Величина вероятности в точности пропорциональна площади пересечения элементарного квадрата и множества решений.

Перейдем теперь к случаю, когда А = (а^) — случайная матрица. Опишем алгоритм на конкретном примере. Для простоты изложения предположим, что ац — случайная дискретная величина, равномерно распределенная на отрезке [2,4]. При этом ац принимает значения 2 + (г — 0.5)Н, г = 1, 2, ...,п; Н = 2/п с вероятностью 1 /п, остальные компоненты совпадают с компонентами матрицы из предыдущего примера. Тогда А — случайная дискретная матрица, принимающая значения Аг с вероятностью рг = 1/п.

Вычислим, к примеру, вероятность Рдх попадания решения в квадрат Дж = [0, 2;0, 3]2. Заметим, что при этом вероятность Рдх складывается из вероятностей попадания решений в Дж для каждой матрицы Аг. Это мы научились делать в предыдущем примере. В данном

случае это будет площадь « фигуры = АгДж П [0,1]2. Окончательно получаем

Рдх = ^ Р®«®.

г

Заметим, поскольку АгДж с [0,1]2, то = |Аг|0.01. После несложных вычислений получаем РдХ = 0,05, что полностью подтверждается численным моделированием методом Монте-Карло. В правой части рис. 5 приведено кусочно-постоянное с шагом 0,1 приближение совместной плотности вероятности вектора ж. Сплошной линией приведена граница множества решений исходной системы.

Рис. 5. Приближение плотности вероятности векторов ж

Рассмотрим пример системы нелинейных уравнений

аж2 + Ьу2 — 4 = 0, жу — с = 0,

где а, Ь, с — равномерные случайные величины, плотности вероятности которых имеют носители [1; 1,1], [2; 2,1], [0, 505; 0, 51]. Зафиксируем с, тогда для каждого ж, у в пространстве параметров 0, 6 подмножество кх у — прямая 0ж2 + 6у2 — 4 = 0. Так несложно вычислить плотность вероятности для вектора решений (ж, у) и далее подсчитать вероятность попадания решения в некоторую область, например, прямоугольник [ж1,ж2] х [у1,у2]. Таким образом, вероятность попадания решения в прямоугольник [0, 37; 0, 375] х [1, 36; 1, 365] — Р = 0,1465. Подобная точность достигается при моделировании методом Монте-Карло с числом бросаний порядка 106.

На рис. 6 приведено кусочно-постоянное приближение совместной плотности вероятности вектора (ж, у) решения системы нелинейных уравнений.

5. Пример использования

В качестве примера использования численных операций над гистограммными переменными рассмотрим задачу принятия решения об инвестировании проекта выпуска лекарственного препарата [9].

Компания рассматривает вопрос о приобретении для последующего производства патента нового лекарственного препарата. Стоимость патента составляет $3,4 млн. Решение

Рис. 6. Приближение плотности вероятности решения системы нелинейных уравнений

принимается на основе анализа дисконтированных денежных потоков NPV и IRR. Горизонт расчетов составляет три года. Стандартная финансовая модель приводится в таблице. Согласно прогнозам, компания в первый, второй и третий годы проекта продаст соответственно 802 тыс., 967 тыс. и 1132 тыс. упаковок лекарства по цене $6, $6,05 и $6,10 за упаковку.

Таблица 1. Стандартная финансовая модель

Год 0 Год 1 Год 2 Год 3

Цена упаковки, $ 6,00 6,05 6,10

Количество проданных шт. 802000 967000 1132000

Выручка, $ 4 812 000 5 850 350 6 905 200

Себестоимость, $, 2 646 600 3 217 693 3 797 860

Валовая прибыль, $ 2 165 400 2 632 658 3 107 340

Операционные издержки, $ 324 810 394 899 466 101

Чистый доход до налогов, $ 1 840 590 2 237 759 2 641 239

Налоги, $ 588 989 716 083 845 196

Стартовые инвестиции, $ -3 400 000

Чистый доход, $ -3 400 000 1 251 601 1 521 676 1 796 043

Ставка налога на прибыль равна 32 %, ставка дисконтирования — 10 %, себестоимость составляет 55%, а операционные издержки — 15% от цены препарата. По результатам расчетов IRR проекта составляет 15 %, а NPV — $344,8 тыс.

В данном случае мы имеем дело с высоким уровнем рыночной неопределенности, поэтому стандартная финансовая модель не может дать достаточных для принятия решения результатов. Для одновременного учета неопределенности в цене, продажах, себестоимости и издержках применяется численный вероятностный анализ. Основные параметры финансовой модели — цена, объем продаж — моделируются как случайные переменные, имеющие вероятностное распределение. Численный вероятностный анализ позволит понять, какие факторы в наибольшей степени повлияют на финансовые результаты проекта. Для моделирования цены продажи (в первый, второй и третий годы проекта отдельно) используется треугольное распределение. Треугольное распределение имеет три параметра — минимальное значение, максимальное значение и наиболее вероятное значение. Цена продажи в первый год имеет минимальное значение $5,90, максимальное значение — $6,10 и наиболее

вероятное значение — $6,00. Аналогично, цена продажи во второй год имеет треугольное распределение с параметрами $5,95; $6,05; $6,15. Цена продажи на третий год имеет треугольное распределение с параметрами $6,00; $6,10; $6,20. Объем продаж моделируется как случайная переменная с нормальным распределением. Объем продаж в первый год имеет нормальное распределение со средним значением (математическим ожиданием) $802 тыс. и стандартным отклонением $25 тыс. Аналогично, объем продаж во второй год имеет нормальное распределение с ожиданием $967 тыс. и стандартным отклонением $30 тыс. Наконец, объем продаж в третий год имеет нормальное распределение с ожиданием $1132 тыс. и стандартным отклонением $25 тыс. Себестоимость (процент от продаж), как предполагается, имеет треугольное распределение с минимальным значением 50 %, максимальным значением 65 % и наиболее вероятным значением 55 %. Следует отметить, что в данном случае треугольное распределение имеет не симметричную форму, а немного скошено вправо, т. е. существует большая вероятность того, что себестоимость будет завышена, а не занижена по сравнению с наиболее вероятным значением. Операционные издержки (процент от продаж) моделируются как нормальное распределение с ожиданием 15 % и стандартным отклонением 2 %.

На рис. 7 приведены гистограммы NPV и IRR для проекта. Из анализа гистограмм NPV и IRR видно, что вероятны как крайне негативные исходы, так и значительная прибыль в сравнении со стандартным анализом. Приведенный пример показывает, что применение ги-стограммной арифметики в рамках технологии визуально-интерактивного моделирования (ВИМ) [13] позволяет ЛПР увидеть возможные варианты негативных исходов реализации проекта, по сравнению со стандартным анализом, который дает только положительный ответ.

-л 1.......rfllll

1 2 10%

Рис. 7. Гистограммы NPV и IRR

20%

-W

30%

0

В заключение отметим, что проведенные авторами теоретические и практические исследования позволяют сделать два основных вывода:

• гистограммная арифметика может рассматриваться как численный метод вероятностного анализа, позволяющий работать с неопределенными данными в рамках различных практических приложений [14];

• гистограммная арифметика может использоваться, как инструмент ВИМ-технологии, что значительно повышает качество анализа возможных вариантов решений и дает в руки ЛПР удобное средство для принятия решений.

Работа выполнена при поддержке гранта ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009-2013 гг., ГК № 02.740.11.0621.

Список литературы

[1] Б.С.Добронец, Интервальная математика, Красноярск, КГУ, 2004.

[2] Л.Т.Ащепков, Д.В.Давыдов, Универсальные решения интервальных задач оптимизации и управления, Институт прикладной математики ДВО РАН, М., Наука, 2006.

[3] А.И.Орлов, Теория принятия решений, Учебное пособие, М., Экзамен, 2005.

[4] В.А.Перепелица, Ф.Б.Тебуева, Дискретная оптимизация и моделирование в условиях неопределенности данных, М., Академия Естествознания, 2007.

[5] Б.С.Добронец, Приближения множеств решений параметрическими множествами, Журнал Сибирского федерального университета. Математика и физика, 2(2009), №3, 305-311.

[6] D.Berleant, Automatically verified reasoning with both intervals and probability density functions, Interval Computations, (1993), №2, 48-70.

[7] W.Li, J.Hym, Computer arithmetic for probability distribution variables, Reliability Engineering and System Safety, 85(2004).

[8] В.А.Герасимов, Б.С.Добронец, М.Ю.Шустров, Численные операции гистограммной арифметики и их применения, АиТ, (1991), №2, 83-88.

[9] А.В.Лукашов, Метод Монте-Карло для финансовых аналитиков: краткий путеводитель, Управление корпоративными финансами, 1(2007), №19, 22-39.

[10] И.М.Соболь, Численные методы Монте-Карло, М., Наука, 1973.

[11] Е.С.Вентцель, Теория вероятностей, М., Наука, 1969.

[12] Б.В.Гнеденко, Курс теории вероятностей, М., Наука, 1988.

[13] Б.С.Добронец, О.А.Попова, Гистограммая арифметика для визуально-интерактивного моделирования в задачах принятия экономических решений, Актуальные проблемы анализа и построения информационных систем и процессов: сб. статей Международной научно-технической конференции, Таганрог, Издательство Технологичсекого института ЮФУ, 2010, 53-57.

[14] Б.С.Добронец, О.А.Попова, Применение гистограммной математики в задачах принятия экономических решении, Труды IX международной ФАМЭТ-2010 конференции, Красноярск, КГТЭН, СФУ, 2010, 127-130.

Numerical Operations on Random Variables and their Application

Boris S. Dobronets Olga A. Popova

It is considers the application of numerical operations on random variables. The description and properties of the basic operations on probability density function are considered. Algorithms and examples of solutions of systems linear algebraic and nonlinear system equations are submitted. It is considered the examples of using of developed methods in practice of economic decision-making. It is shown, that the given approach allows essentially improve quality of decision-making and reduce volume of calculations.

Keywords: numerical operations on random variables, histogram arithmetics, functions of random variables, stochastic system linear algebraic equations, stochastic nonlinear equations, the theory of decision-making, interval mathematics.