Моделирование случайной величины методом суперпозиции

[probably]

Дана плотность случайной величины $\xi$: $p(x) = c\ln^4 x, x \in \left( 0,1\right]$. Сначала рассмотрим аналогичную задачу для с.в. $\gamma$: $p(x) = -\ln x, x \in \left( 0,1\right]$. Если сделать замену $y=1-x$, то плотность станет $p(y) = -\ln(1-y), y \in \left[0,1\right) $. Разложив логарифм в ряд, получаем $p(x) = \sum\limits_{k=1}^\infty \frac{y^k}{k} = \sum\limits_{k=1}^\infty \frac{1}{k(k+1)}(k+1)y^k$. Получили смесь распределений $f_{k}(y) = (k+1)y^k$ с весами $p_{k} = \frac{1}{k(k+1)}$. Теперь, чтобы смоделировать $\gamma$, нужно сначала смоделировать вспомогательную дискретную с.в. $\eta$ со значениями $p_{k}$, взять получившееся $k$, подставить в $f_{k}(y)$ и оттуда найти $\gamma$. Тогда $\int\limits_{0}^{\gamma}(k+1)y^k dx = \alpha$, где $\alpha$ распределено равномерно на [0,1]. Т.е. $y^{k+1}\rvert\limits_{0}^{\gamma} = \alpha$, или $\gamma = \alpha^\frac{1}{k+1}$. Теперь вернемся к исходной задаче. Сначала найдем константу $c$ из условия, что интеграл от плотности равен единице. $c = \frac{1}{24}$. Сделаем замену $y = 1 - x$. Отделим один логарифм и разложим его в ряд: $p(x) = \sum\limits_{k=1}^\infty \frac{1}{k(k+1)}(k+1)y^{k}\ln^3(1-y)$. Тогда задача сведется к моделированию с.в. $\xi_{1}$ с плотностью $p_{1}(x) = (k+1)y^k\ln^3(1-y)$, где $k$ уже фиксирована. Аналогично вышеописанному раз за разом раскладываем логарифмы в ряд. В итоге получим, что надо реализовать случайную величину $\xi_{4}$ с плотностью $p_{4}(x) = \frac{1}{24}(k+1)y^k(m+1)y^m(n+1)y^n(s+1)y^s$ или $p_{4}(x) = \frac{1}{24}(k+1)(m+1)(n+1)(s+1)y^{k+m+n+s}$, где значения $k, n, m, s$ известны. Т.е. $\xi = \left(\frac{24(k+m+n+s+1)\alpha}{(k+1)(m+1)(n+1)(s+1)}\right)^{\frac{1}{k+m+n+s+1}}$. Проблема в том, что у полученной случайной величины мат. ожидание отличается от теоретического знаком. Думаю, что дело в замене, но только для $p(x) = -\ln x$ получается всё правильно (мат. ожидание и дисперсия выходят с правильным знаком). Не могу сообразить, в чём косяк...

[--mS--]

probably в сообщении #662505 писал(а):

с плотностью $p_{1}(x) = (k+1)y^k\ln^3(1-y)$,

Логарифм отрицателен. Минус напишите.

А по поводу матожиданий непонятно. Матожидание получается отрицательным? Как бы Вы ни считали матожидание, это невозможно: все величины, что Вы моделируете, принимают значения на $(0,\,1)$.

Что же до замены, то она должна катастрофически влиять и на случай одного логарифма. Смотрите, что Вы делаете: плотность $p_\xi(1-x)$ есть производная от функции $1-F_\xi(1-x)$, т.е. это плотность случайной величины с функцией распределения

$$1-F_\xi(1-x) = 1-\mathsf P(\xi < 1-x)=\mathsf P(\xi > 1-x) = \mathsf P(1-\xi < x).$$

Т.е. Вы моделируете не $\xi$ с плотностью $c\ln^4 x$ или $-\ln x$, а величину $1-\xi$. Соответственно, для логарифма в первой степени матожидание у моделируемой величины будет $3/4$ вместо $1/4$, как должно быть у $\xi$. То же самое будет и с четвертой степенью логарифма: вместо матожидания $1/32$ получится $31/32$.

[probably]

С минусом опечатка, у нас он есть.

Мат.ожидание получается отрицательным потому, что у нас генерируются почему-то и отрицательные случайные велечины (где-то ошибка в логике видимо, которую тоже надо найти).

Замена понятно что влияет. Забыл написать, что мы в итоге берем случайную величину: 1 минус наша сгенерированная величина (т.е. в конечной формуле для кси забыли написать "1-...").

Так что всё равно не понятно в чём именно проблема.

[--mS--]

Это как это арифметический корень (который генерируется) может оказаться отрицательной с.в.???

probably в сообщении #662505 писал(а):

$\xi = \left(\frac{24(k+m+n+s+1)\alpha}{(k+1)(m+1)(n+1)(s+1)}\right)^{\frac{1}{k+m+n+s+1}}$.

[aaaaa]

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

[probably]

Правильно ли я рассуждаю:
Если положить $\eta = 1 - \xi$, то $p_{\eta}(x) = p_{\xi}(1-x)$, и для нахождения реализации $\eta$ нужно считать $\alpha = \int\limits_{0}^{\eta}p_{\eta}(y)dy = \int\limits_{1}^{1-\eta}p_{\xi}(1-x)d(1-x) = \int\limits_{1-\eta}^{1}p_{\xi}(1-x)dx$? И можно ли продолжить равенство дальше $ = \int\limits_{\xi}^{1}p_{\xi}(1-x)dx$, так что сразу находится реализация $\xi$?

[--mS--]

Правильно, конечно. Вот теперь наконец понятно, раз больше единицы у Вас получаются интегралы от плотности. Проблема вот где: после отщепления одного логарифма получается ряд с коэффициентами, сумма которых единица, вот только сами члены этого ряда $-\frac{k+1}{24}x^k\ln^3 x$ плотностями не являются. Плотностью будет $-\frac{(k+1)^4}{6}x^k\ln^3 x$. К сожалению, ряд из коэффициентов, которые вылезут, если выделять такие плотности, а именно, ряд $\sum_{k=1}^\infty \frac{1}{4k(k+1)^4}$, единице не равен.

Ну и соответственно, плотность, по которой Вы получаете итоговую $\eta$ (или $\xi$) - которая корень страшной степени из альфы с иногда большим множителем, - это совсем не плотность, вот поэтому её интеграл оказывается больше единицы, что потом и даёт минус.

Что с этим делать и можно ли как-то исправить, не знаю.

(Оффтоп)

Однозначно очень большое спасибо за подробное и понятное изложение метода. Я как-то далека от методов моделирования с.в., но с удовольствием вникла, Вы хорошо пишете.

[probably]

Интересно, что такой способ моделирования даёт почти правильный ответ при моделировании плотности $p(x) = \frac{1}{2}\ln^{2}(x)$ (мат. ожидание отличается примерно на 0.12, а дисперсия получается такая же), при моделировании $p(x) = \frac{1}{24}\ln^{4}(x)$ мат. ожидание и дисперсия отличаются лишь знаками (по модулю - незначительно), а для $p(x) = \frac{1}{720}\ln^{6}(x)$ мат. ожидание и дисперсия отличаются от теоретической весьма значительно... Повезло, наверное.

[--mS--]

Ну вообще функция распределения с плотностью $c\ln^4x$ выписывается, можно численно её обращать в нужной точке $\alpha$, являющейся реализацией равномерно распределённой с.в., так что в принципе моделировать нужную с.в. нет проблем. Но вот как и можно ли поправить такой приятный метод (люблю линейные комбинации распределений) - любопытно.