CC BY
133
Естественные науки
УДК 517.518.8
АППРОКСИМАЦИЯ ПЛОТНОСТЕЙ РАСПРЕДЕЛЕНИЙ СЛУЧАЙНЫХ ВЕЛИЧИН С ПРИМЕНЕНИЕМ ОРТОГОНАЛЬНЫХ ПОЛИНОМОВ ЧЕБЫШЕВА-ЭРМИТА
А.В. Димаки, А.А. Светлаков
Томский государственный университет систем управления и радиоэлектроники E-mail: [email protected]
Предложен способ аппроксимации плотностей статистических распределений полиномами Чебышева-Эрмита для дальнейшей генерации последовательностей случайных чисел, закон распределения которых близок к аппроксимируемому закону. Разработана методика регуляризации решения системы линейных алгебраических уравнений, возникающей при решении задачи аппроксимации плотности распределения. Исследованы возможности применения предложенного способа для аппроксимации гладких и островершинных распределений, имеющих одну моду, а также двухмодальных распределений.
1. Введение
При моделировании поведения какого-либо процесса или объекта, находящегося под воздействием случайных факторов, возникает задача получения случайных чисел, закон распределения которых близок к закону распределения, наблюдаемому в эксперименте [1]. В этом случае зачастую не важно знать тип и значения параметров закона распределения генерируемых случайных чисел, важно лишь, чтобы этот закон соответствовал закону распределения, полученному при наблюдении за поведением исследуемого объекта или процесса.
Один из подходов к решению указанной задачи заключается в нахождении некоторой аппроксимирующей функции, которая описывает экспериментальный закон распределения с требуемой точностью. На основе этой функции искомые случайные числа могут быть получены из равномерных случайных чисел методом обратной функции [2].
Несмотря на кажущуюся простоту такого подхода, при его практической реализации возникает ряд проблем, связанных, в первую очередь, с тре-
бованиями, которые предъявляются к функциям, аппроксимирующим плотность распределения и функцию распределения. К этим требованиям относятся: 1) положительность функции, аппроксимирующей плотность распределения, на всей ее области определения и 2) равенство интеграла, взя-
того от аппроксимирующей функции по всей области ее определения, единице. Кроме того, зачастую возникают проблемы, связанные с низкой устойчивостью получаемых решений. Одним из способов, позволяющих успешно решить задачу аппроксимации плотности распределения, является подход, основанный на применении ортогональных полиномов Чебышева-Эрмита [3].
2. Постановка и решение задачи аппроксимации плотности распределения полиномами Чебышева-Эрмита
Как известно, многочлены Чебышева-Эрмита Н(х) ортогональны на всей числовой оси х с весовой функцией
h( x) = e
(1)
и, соответственно, удовлетворяют соотношениям
» (О,I ф у -
I" Н (х)Н (х)Н(х)с1х = \ , I,у = 0...<».
I5,' = у
В случае, когда используются ортонормирован-ные многочлены Чебышева-Эрмита, значение переменной я равно единице.
Вычисление многочленов Чебышева-Эрмита может осуществляться при помощи формулы Ро-дрига, имеющей вид [3]:
х
Ип(х) = (-1)пех (е)(п>,
либо при помощи рекуррентного соотношения вида Ии+1(х) = 2хИп (х) - 2иИи-1(х).
Как показано в работе [3], частичные суммы рядов Фурье для некоторой функции у(х) по ортогональным многочленам являются наилучшими приближениями данной функции на конечном интервале [а; Ь]. Тот факт, что весовая функция для ортогональных многочленов Чебышева-Эрмита является решением известного уравнения Пирсона [3], делает предпочтительным использование данного семейства многочленов для аппроксимации плотностей распределений с учетом того, что они ортогональны на всей числовой оси. Следует отметить также, что при разложении функции у(х) в ряд Фурье по многочленам Чебышева-Эрмита к ней предъявляются минимальные требования в смысле ее интегрируемости.
В рамках данной работы аппроксимация плотности распределения, наблюдаемой в эксперименте, производилась при помощи функции следующего вида:
У (х) = Ё а,И, (х),
(2)
* = Ё
]=1
у] -Ё а,И, (х)
^ Ш1П,
(3)
где т - число интервалов гистограммы, построенной по экспериментальной выборке.
Дифференцируя выражение (3) по неизвестным коэффициентам ак, получаем систему из п линейных уравнений следующего вида:
да,.
= -Ё
]=1
у. -Ё а,И, (х)
Ик (X.), к = 0...п. (4)
После ряда очевидных преобразований система уравнений (4) может быть представлена в следующей матричной форме:
Ва = у, (5)
где а - вектор искомых коэффициентов а,, а элементы Ьц матрицы В и компоненты у1 вектора у вычисляются следующим образом:
Ь . =Ё И (хк И (хк),
к=1
т
У> = Ё укИ,(хк)-
(6)
(7)
Очевидно, что матрица В является симметричной, и ее элементы, расположенные на главной диагонали, отличны от нуля. Следует отметить также, что элементы матрицы В, не расположенные на
главной диагонали, также в общем случае отличны от нуля, что обусловлено следующими факторами:
1. в выражении (6) при суммировании полиномов не используется весовая функция (1);
2. интервал значений переменной х, на котором производится суммирование, конечен;
3. суммирование производится только в конечном множестве точек хк, к=1,...,т.
Решая систему (5), можно найти вектор искомых коэффициентов а. Предложенный подход позволяет с необходимой точностью аппроксимировать плотности одномодальных распределений, форма которых близка к форме нормального распределения (рис. 1).
где п - максимальный порядок многочлена; а1 -постоянные коэффициенты.
В этом случае задача аппроксимации сводится к нахождению таких значений коэффициентов ак, при которых выполняется условие п ~\2
Рис. 1. Аппроксимация плотности р(х) нормального распределения при помощи полиномов Чебышева-Эрмита
Однако осуществить аппроксимацию плотностей распределений, форма которых значительно отличается от нормального распределения, в частности двухмодальных распределений, не удается. В частности, наблюдается появление «всплесков» аппроксимирующей функции. Кроме того, на краях области определения в некоторых случаях возникают участки, где аппроксимирующая функция опускается ниже оси абсцисс. Примером такого поведения описанного алгоритма являются результаты аппроксимации плотности двухмодального распределения (рис. 2).
Как показано на рис. 2, варьирование максимального порядка аппроксимирующего полинома не позволяет осуществить аппроксимацию таким образом, чтобы выполнялись сформулированные выше требования. Данный факт обусловлен, по-видимому, слабой устойчивостью решения системы (5) к ошибкам счета, возникающих при вычислении полиномов, входящих в сумму (2), причем, очевидно, между порядком этих полиномов и величинами ошибок имеется некоторая взаимосвязь.
Для повышения устойчивости решения системы уравнений (5) авторами был применен широко известный метод регуляризации, предложенный А.Н. Тихоновым [4]. В рамках данного метода точное решение а системы уравнений (5) заменяется некоторым приближенным решением другой системы, которая имеет следующий вид:
(В + /Е) а = у, у> 0, (8)
где Е - единичная матрица, у-некоторый малый коэффициент, называемый параметром регуляризации.
1=0
к=1
Рис. 2. Результаты аппроксимации плотности р(х) двухмодального распределения при различных значениях порядка аппроксимирующего полинома п
Рассматривая эту систему, можно непосредственно видеть, что при: 1) больших значениях у решение ее сильно отличается от решения системы (5), однако является устойчивым; 2) малых уреше-ние близко к решению системы (5), однако устойчивость его уменьшается. Следовательно, необходимо выбирать такое значение у, при котором точность и устойчивость получаемого решения являются оптимальными [4].
Как показали проведенные исследования, в общем случае невозможно подобрать такое значение у, при котором обеспечивалась бы достаточная точность аппроксимации при отсутствии колебаний аппроксимирующей кривой на краях ее области определения. Данное явление связано, очевидно, с существенно различными значениями ошибок счета, возникающих при вычислении полиномов различных степеней, и, как следствие, с различными значениями ошибок задания элементов матрицы В и вектора у в системе (5). Таким образом, представляется логичным заменить в (8) единичную матрицу Е на некоторую матрицу весовых коэффициентов М, сформированную на основании имеющейся информации о величинах ошибок вычисления каждого из полиномов, либо хотя бы о соотношении значений этих ошибок. При этом система (8) преобразуется к следующему виду:
(В + аМ)а = у, а > 0, (9)
где а - малый положительный параметр.
Известно, что ошибки вычислений на ЭВМ, проводимых с использованием чисел с плавающей точкой, определяются множеством факторов, влияние которых в большинстве практических случаев учесть невозможно [5]. В качестве примеров подобных факторов можно привести особенности архитектуры центрального процессора, используемой
операционной системы, языка программирования, версии компилятора и т. д. В связи с этим одним из возможных способов формирования матрицы М, входящей в систему уравнений (9), является анализ значений элементов матрицы В и вектора у.
Пусть каждый из полиномов Щ(х) вычисляется в точке хк с ошибкой АЦ. Тогда выражения (6, 7) можно переписать следующим образом:
т
Ъи = £ (Н (Хк ) + Щ )(И] (хк ) + М/ ) =
к=1 т (
= £
Н, (Хк ) Н] (Хк ) + Н. (хк )Щ +
(хк )АЙ/ + Щ АН'к
т
У = £ Ук (Н (хк) + К).
(10)
(11)
Разделив обе части равенства (10) на соответствующие части равенства (11), пренебрегая в выражении (10) членом высшего порядка малости, получаем:
Ъ]
У
£ (Н х )Н] х) + Н] х )Щ + Н (хк )АК])
. к=1_
т
£ ^к (Н (хк)+К)
(12)
При 1=] выражение (12) преобразуется к следующему виду:
£ (Н (хк )(Н (хк) + 2АК)) £Н\х) Ъ± = _~ _(13)
У т т
£ ук (Н (хк) + К) £ УН (хк)
к=1
к= 1
Важными особенностями отношения (13) являются следующие:
к=1
к=1
к=1
1. при стремлении ошибки Ак| к нулю отношение (13) стремится к отношению точных значений элементов Ьц и у;
2. в зависимости от знака ошибки отношение (13) может отличаться от отношения истинных значений элементов Ь ц и у; как в меньшую, так и в большую сторону;
примерное равенство в (13) справедливо только при I = ], т. к. значения ошибок вычисления полиномов различных степеней могут сильно различаться.
На основании сказанного выше, сформируем матрицу М следующим образом:
Ъ,
3
у, увеличение параметра а не столь негативно сказывается на точности аппроксимации. При а>1 наблюдается подавление колебаний аппроксимирующей кривой, и вместе с тем точность аппроксимации снижается незначительно. В целом, наилучшие результаты аппроксимации в смысле максимальной точности и минимальных колебаний были получены при значениях ае (2;5).
Важным достоинством предложенного алгоритма является высокая устойчивость получаемых решений при увеличении порядка п полинома (2). На рис. 4 приведены аппроксимирующие кривые, полученные при п=50.
М = {т] =
' = ]
У1 0,,* ]
Определим теперь значение малого коэффициента а, входящего в выражение (7). Для удобства его определения пронормируем матрицу М следующим образом:
М = {т] =
т,
' = ]
тах(т,,)' 0,, * ]
Таким образом, значения элементов матрицы М будут находиться в интервале от нуля до единицы.
Рассмотрим влияние значения параметра а на результаты аппроксимации при различных порядках полинома п (рис. 3).
Как видно из рис. 3, варьируя значение а, можно добиться достаточной точности аппроксимации при отсутствии «всплесков» и отрицательных значений аппроксимирующей кривой. Необходимо отметить, что, по сравнению с влиянием параметра
Рис. 4. Результаты аппроксимации плотности p(x) двухмо-дального распределения при различных значениях параметра а и п = 50
Как показывают результаты, полученные при аппроксимации экспериментальной двухмодальной плотности распределения, при увеличении п точность аппроксимации увеличивается (рис. 5). Без использования регуляризации точность аппроксимации максимальна, однако с ростом п наблюдаются существенные колебания аппроксимирующей функции на краях области ее определения. Применение метода регуляризации А.Н. Тихонова [4] позволяет подавить колебания аппроксимирующей функции, но точность аппроксимации при этом
снижается в 5... 6 раз. Применение изложенного в данной работе алгоритма регуляризации позволяет подавить колебания аппроксимирующей функции и получить точность аппроксимации, сравнимую с той, которая имеет место без использования регуляризации. Таким образом, результаты, приведенные на рис. 4 и 5, являются подтверждением высокой устойчивости решений, получаемых с помощью разработанного алгоритма, к ошибкам вычисления полиномов Чебышева-Эрмита высоких степеней.
п
Рис. 5. Зависимости среднеквадратического отклонения 5 от порядка полинома п, построенные при аппроксимации плотности распределения, наблюдаемого в эксперименте, с использованием регуляризации и без нее
3. Применение результатов аппроксимации плотностей распределений для получения последовательностей случайных чисел с заданным законом распределения
Замечательной особенностью разработанного способа аппроксимации плотностей распределений, основанного на применении полиномов Чебышева-
Эрмита, является возможность аналитического вычисления интеграла от аппроксимирующей функции. В работе [3] показано, что определенный интеграл от полинома Чебышева-Эрмита к-го порядка Нк(х) может быть представлен следующим образом:
к < **=■
а
Следовательно, определенный интеграл от аппроксимирующей функции (2), взятый на интервале [х1; х2], может быть вычислен при помощи простого выражения:
] у(Х)ёХ = £ а, Н+'(^• (14)
Таким образом, для вычисления интеграла от аппроксимирующей функции не требуется применения каких-либо численных методов. Отметим, что, поскольку функция (2) аппроксимирует плотность некоторого распределения, то интеграл от (2), вычисленный при помощи выражения (14), аппроксимирует, очевидно, функцию данного распределения. Эта особенность позволяет при генерации случайных чисел, закон распределения которых близок к экспериментальному закону распределения, использовать высокопроизводительные методы вычисления обратной функции, такие как метод дихотомии, метод Ньютона и др.
На рис. 6 приведены гистограммы, аппроксимирующие их функции, а также интегралы от аппроксимирующих функций, построенные для модельных выборок и выборок значений случайных величин, наблюдаемых в эксперименте.
Как видно из рис. 6, удовлетворительные результаты аппроксимации для различных законов распределения получаются при существенно раз-
-2-10 1 X
а) и = 8; 5"= 0,003
-1 -0,5 0 0,5 1
в) п = 35; 5 = 0,010
б) п = 10; 5 =0,021
-2-10 1 2 г) п= 33; 5 = 0,020
Рис. 6. Гистограммы, аппроксимирующие их функции, а также интегралыы от этих функций; а) нормальное распределение; б) гамма-распределение; в) бета-распределение; г) распределение Коши; а=2
личающихся порядках полинома п. Кроме того, аппроксимация островершинных распределений, таких как распределение Коши, не может быть осуществлена с достаточно высокой точностью. Увеличение п в подобных случаях не приводит к улучшению качества аппроксимации, т. к. при этом увеличиваются колебания аппроксимирующей функции. Результаты, приведенные на рис. 6, позволяют сделать вывод о том, что область применения разработанного подхода ограничена; наилучшим образом данный подход может быть использован при аппроксимации достаточно гладких плотностей распределений.
Несмотря на этот факт, интеграл от аппроксимирующей функции (2) близок к интегралу от аппроксимируемой плотности распределения. Таким образом, возможно применение разработанного подхода при генерации случайных чисел, подчи-
няющихся требуемому закону распределения, методом обратной функции.
Заключение
Разработанный способ аппроксимации плотностей распределений на основе ортогональных полиномов Чебышева-Эрмита может быть с успехом применен при обработке экспериментальных выборок и перспективен при решении широкого круга научных и практических задач, связанных с моделированием случайных величин. Предложенная методика регуляризации позволила получить устойчивые решения даже при достаточно высоких порядках (п<50) аппроксимирующих полиномов. Полученная аппроксимирующая функция может быть использована при создании эффективных алгоритмов генерации случайных чисел с заданным законом распределения.
СПИСОК ЛИТЕРАТУРЫ
1. Ермаков С.М. Метод Монте-Карло и смежные вопросы. Изд. 2-е. - М.: Наука, 1975. - 412 с.
2. Бобнев М.П. Генерирование случайных сигналов. Изд. 2-е. -М.: Энергия, 1971. - 290 с., ил.
3. Суетин П.К. Классические ортогональные многочлены. - М.: Наука, 1976. - 328 с.
4. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. - М.: Наука, 1974. - 186 с.
5. Мудров А.Е. Численные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль. - Томск: РАСКО, 1992. - 272 с.
УДК 519.22(075.8)
СИНТЕЗ И НЕКОТОРЫЕ РЕЗУЛЬТАТЫ ИССЛЕДОВАНИЙ АЛГОРИТМОВ РЕШЕНИЯ УРАВНЕНИЙ, ВОЗНИКАЮЩИХ ПРИ ИСПОЛЬЗОВАНИИ РАСПРЕДЕЛЕНИЯ СТЬЮДЕНТА
А.А. Светлаков, Ю.Г. Свинолупов*, Е.В. Шумаков
Томский государственный университет систем управления и радиоэлектроники *ОАО «Манотомь», г. Томск E-mail: [email protected]
Исследованы численные методы решения уравнений, возникающих при использовании t-распределения для интервального оценивания характеристик случайных величин. Приведены результаты решения уравнения с помощью метода Ньютона, хорд и, наиболее оптимального, дихотомического деления.
1. Введение
В настоящее время уделяется большое внимание разработке автоматизированных систем управления технологическими процессами (АСУТП). Одной из главных задач, возникающих в процессе функционирования подобных систем, является задача получения достаточно точных оценок неизвестных истинных характеристик управляемого объекта [1]. Во многих практических случаях при определении данных характеристик объектов для обеспечения желаемой точности целесообразно определять не точечные оценки их неизвестных значений, а некоторый интервал, который с заданной вероятностью накрывает данные значения. Обычно при решении задач по определению подобных оценок, особый интерес представляет не только точность, но и надежность
оценок при ограниченном числе измерений. В математической статистике решение этих задач производится путем построения доверительного интервала 1а, который с заданной вероятностью а накрывает оцениваемое значения [1-4].
Доверительный интервал 1а устанавливает точность определения оценок рассматриваемого параметра, а характеристикой надежности получаемой оценки является доверительная вероятность а. Существующий способ построения доверительных интервалов с помощью заранее составленных и широко известных в математической статистике таблиц малопригоден для применения в АСУТП, во-первых, из-за ограниченности значений доверительной вероятности, для которых составлены данные таблицы, и во-вторых, их использование связа-
