Генератор многомерного распределения Коши

[Sinclair]

Добрый день.

Мне нужно написать генератор многомерного распределения Коши, функция плотности вероятности для которого выглядит так (нормировку опускаю, центр распределения - начало координат):
$\text{PDF}(\vec{x})=\frac{T}{\left(\vec{x}^2+T^2\right)^{\frac{D+1}{2}}}$,
где:
$D \in \mathbb N$ - число измерений.
$T \in \mathbb R, T > 0$ - параметер масштабирования (scale).
$\vec{x} \in \mathbb R^D$ - радиус-вектор в $D$-мерном пространстве.
$\vec{x}^2$ - скалярное произведение $\vec{x}$ на самого себя.
Я не смог нигде найти генератор для многомерной версии этого распределения, так что, мне пришлось написать свой.
К сожалению, я неуверенно чувствую себя в этой области, а потому, хочу проверить свое решение:
Поскольку, компоненты $\vec{x}$ зависимы, то мы не можем сгенерировать вектор покомпонентно.
Однако, можно заметить, что для всех векторов одинаковой длины (то есть, векторов, у которых совпадает $\vec{x}^2$), функция дает одинаковую плотность вероятности. Таким образом, в полярных координатах, плотность вероятности зависит только от длины $\vec{x}$, и никак не зависит от направления.
Это позволяет нам разбить генерацию вектора на 2 части:
1) Генерацию направления вектора, что тривиально сводится к генерации равномерно распределенных углов для полярных координат, и не представляет дальнейшего интереса.
2) Генерацию длины вектора.
Генерацию длины вектора я полагаю осуществить следующим образом:
2a) Получить фукнцию плотности вероятности для длины вектора.
2b) Получить функцию распределения вероятности, проинтегрировав плотность вероятности.
2c) Найти обратную функцию к функции распределения вероятности.
Тогда, генерировать длину вектора можно будет просто выбрав равномерно распределенную случайную величину из области значений функции распределения вероятности, и, получив длину с помощью обратной функции.
Теперь подробнее:
2a) Нужно перейти от плотности вероятности для вектора, к плотности вероятности для длины вектора.
Как мне подсказывает математическая интуиция, и соображения о пределе отношения объема полых гиперсфер с большим и малым радиусом, плотность вероятности длины вектора получается домножением исходной плотности вероятности в произвольной точке с данной длиной на площадь поверхности гиперсферы с радиусом, равным этой длине (иначе говоря - на интеграл по всей области с данной длиной).
Языком формул, получается следующий результат (поскольку, мне не нужен нормировочный коэффициент, от радиуса сферы остается только $r^{D-1}$):
$\text{PDF}_r(r)=\frac{r^{D-1}T}{\left(r^2+T^2\right)^{\frac{D+1}{2}}}$
К сожалению, мои математические познания не достаточны для того, чтобы твердо доказать этот способ, и я даже не знаю, как это можно проверить.
2b) Теперь, имея плотность вероятности, я могу найти функцию распределения вероятности длины вектора.
Для этого мне нужно взять интеграл $\sideset{}{_0^Y}\int\text{PDF}_r(r)\text{dr}$. Сам взять подобный интеграл я вряд ли смогу, я взял его с помощью Wolfram Mathematica. Результат:
$\frac{\left(\frac{Y}{T}\right)^D\sideset{_2}{_1}F\left[\frac{D}{2},\frac{1+D}{2},\frac{2+D}{2},-\frac{Y^2}{T^2}\right]}{D}$,
где $\sideset{_2}{_1}F$ - гипергеометрическая функция. С помощью, опять же, вольфрама, я могу найти предел этой функции при $Y$, стремящимся к бесконечности, для определеня границ генерации равномерно распределенного случайного значения.
2c) После этого, остается найти обратную функцию, и генератор готов.
Здесь вознакают серьезные трудности, потому что нужно найти обратную фукнкцию к произведению показательной и гипергеометрической функций. Как это сделать я не представляю. К счастью, этого можно избежать. Первый способ - построить обратную фукнцию с помощью кусочно-линейной аппроксимации. Второй - используя то, что функция распределения монотонно возрастает, каждый раз искать обратное значение с помощью двоичного поиска.
Отсюда возникает ряд вопросов.
Первый и самый главный - правильно ли я вообще все делаю?
То есть, нет ли у меня структурной ошибки в написании генератора?
Если вы разбираетесь в теме, и видите, что все правильно, пожалуйста, напишите это явно.
Второй - можно ли, и если да, то как, найти в общем виде обратную функцию на шаге 2c? Не обязательно точное значение, подойдет также ряд, или любой способ получить решение с задаваемой точностью. Очевидно, что это способ должен быть быстрее двоичного поиска, иначе в нем нет смысла.
Третий. Я еще сам в этом направлении не копал, но задам вопрос, раз уж пишу - какие вы могли бы порекомендовать библиотеки для вычисления гипергеометрической функции? (ответ имеет смысл только если нельзя найти обратную функцию в общем виде без гипергеометрической). Пишу на Cпп под линукс, главное - скорость.

[--mS--]

До п. 2а включительно всё верно. Плотность распределения длины вектора именно такая. Дальше - не знаток.

Думаю, однако, что есть способ генерации проще. Ваше многомерное распределение Коши есть частный случай многомерного же распределения Стьюдента с одной степенью свободы: с диагональной матрицей $\Sigma=T^2 \mathbb E$.

Пусть $Y_1,\ldots,Y_D,\,Y_{D+1}$ - независимые стандартные нормальные случайные величины. Тогда вектор $\vec Y = (TY_1,\ldots, TY_D)$ имеет нормальное многомерное распределение с нулевым вектором средних и матрицей ковариаций $\Sigma=T^2\mathbb E$, и отношение $\displaystyle \frac{\vec Y}{Y_{D+1}}$ имеет как раз искомое многомерное распределение Коши.

[Sinclair]

--mS-- в сообщении #588409 писал(а):

Пусть $Y_1,\ldots,Y_D,\,Y_{D+1}$ - независимые стандартные нормальные случайные величины. Тогда вектор $\vec Y = (TY_1,\ldots, TY_D)$ имеет нормальное многомерное распределение с нулевым вектором средних и матрицей ковариаций $\Sigma=T^2\mathbb E$, и отношение $\displaystyle \frac{\vec Y}{Y_{D+1}}$ имеет как раз искомое многомерное распределение Коши.

Да, это действительно ответ на мой вопрос.
Спасибо.
Значит вчера занимался ерундой.
Еще один день на свалку.

[--mS--]

Да ну почему ерундой, всё равно полезно то, что Вы проделали: не здесь, так в ином месте пригодится.

[Sinclair]

--mS-- в сообщении #588423 писал(а):

Да ну почему ерундой, всё равно полезно то, что Вы проделали: не здесь, так в ином месте пригодится.

Это дежурное утешение для учащихся, чтобы мотивировать их не сдаваться.
На самом деле, это не так. Ярчайший пример - 90% школьной программы.

-- 24.06.2012, 11:21 --

--mS-- в сообщении #588409 писал(а):

многомерное распределение Коши есть частный случай многомерного же распределения Стьюдента с одной степенью свободы: с диагональной матрицей $\Sigma=T^2 \mathbb E$.

Пусть $Y_1,\ldots,Y_D,\,Y_{D+1}$ - независимые стандартные нормальные случайные величины. Тогда вектор $\vec Y = (TY_1,\ldots, TY_D)$ имеет нормальное многомерное распределение с нулевым вектором средних и матрицей ковариаций $\Sigma=T^2\mathbb E$, и отношение $\displaystyle \frac{\vec Y}{Y_{D+1}}$ имеет как раз искомое многомерное распределение Коши.

и вот еще один вопрос.
Здесь вы написали, что домножение единичной матрицы ковариации на некоторую константу в квадрате аналогично домножению сгенерированных единичной матрицей величин на эту константу. Насколько я понимаю, домножение матрицы аналогично изменению множителя аргумента экспоненты.
Это ошибка, или это свойство нормального распределения?
Понятно, что дисперсия изменится в соответствии с именно этим законом, но я не могу понять, изменится ли остальное распределение аналогично (то есть, равно ли изменение дисперии растягиванию/сжатию гауссианы).

[--mS--]

Sinclair в сообщении #588425 писал(а):

На самом деле, это не так. Ярчайший пример - 90% школьной программы.

Если 90% школьной программы выкинуть, затем можно смело выкинуть освоивших остальные 10%.

Sinclair в сообщении #588425 писал(а):

Здесь вы написали, что домножение единичной матрицы ковариации на некоторую константу в квадрате аналогично домножению сгенерированных единичной матрицей величин на эту константу. Насколько я понимаю, домножение матрицы аналогично изменению множителя аргумента экспоненты.
Это ошибка, или это свойство нормального распределения?
Понятно, что дисперсия изменится в соответствии с именно этим законом, но я не могу понять, изменится ли остальное распределение аналогично (то есть, равно ли изменение дисперии растягиванию/сжатию гауссианы).

Прошу простить, не поняла вопрос. Нормальное распределение определяется вектором средних и матрицей ковариаций. Так же как и одномерное нормальное определяется матожиданием и дисперсией. Если у независимых стандартных нормальных величин матрица ковариаций единична, то у тех же величин, умноженных на константу, дисперсии возрастут в квадрат константы раз, и матрица ковариаций станет диагональной с этой константой в квадрате на диагонали. Нормальность же распределения при (невырожденных) линейных преобразованиях векторов не портится. Соответственно, $T\vec X \sim \textrm N_{\vec 0, \Sigma}$, где $\Sigma=T^2\mathbb E$. Плотность при этом будет равна

$$f_{T\vec X}(\vec x) = \frac{1}{\sqrt{\det \Sigma}(2\pi)^{D/2}} \exp\left(\frac12 (\vec x - \vec \mu)^T\Sigma^{-1}(\vec x - \vec \mu)\right)= \frac{1}{T^D(2\pi)^{D/2}} \exp\left(-\frac12\sum_{i=1}^D \frac{x_i^2}{T^2}\right).$$

Первое равенство - просто определение плотности нормального распределения с заданной матрицей ковариаций $\Sigma$ и вектором средних $\vec \mu$.

Да в конце концов тут и не нужно никаких многомерностей, компоненты независимы, достаточно выяснить распределение одной компоненты. А если $X\sim \textrm N_{0,1}$, то $TX\sim N_{0, T^2}$ - просто потому, что нормальность распределения не портится, а дисперсия изменилась в $T^2$ раз. Вот и плотность будет $f_{TX}(x)=\frac{1}{T}f_X\left(\frac{x}{T}\right)=\frac{1}{T\sqrt{2\pi}}\exp\left(-\frac{x^2}{2T^2}\right)$.

[Sinclair]

--mS-- в сообщении #588478 писал(а):

Если 90% школьной программы выкинуть, затем можно смело выкинуть освоивших остальные 10%.

Можно выкинуть (освоивших), конечно, но ведь можно и не выкидывать же.

--mS-- в сообщении #588478 писал(а):

Да в конце концов тут и не нужно никаких многомерностей, компоненты независимы, достаточно выяснить распределение одной компоненты. А если $X\sim \textrm N_{0,1}$, то $TX\sim N_{0, T^2}$ - просто потому, что нормальность распределения не портится, а дисперсия изменилась в $T^2$ раз. Вот и плотность будет $f_{TX}(x)=\frac{1}{T}f_X\left(\frac{x}{T}\right)=\frac{1}{T\sqrt{2\pi}}\exp\left(-\frac{x^2}{2T^2}\right)$.

В общем, я это и хотел выяснить - испортится ли нормальность одномерного распределения, при домножении сгенерированных величин на константу.
Я, кстати, так и не понимаю, из чего это следует, что нет.
С другой стороны, если генерировать случайные величины предложеным мною методом - выбором случайной величины с области значений функции распределения (интеграла от плотности), и вычислением обратной - то да, в этом случае домножение аргумента экспоненты на $T^2$ эквивалентно увеличению сгенерированных значений в $T$ раз, поскольку при интегрировании эта константа становится множителем к аргументу функции ошибки.
Но это все интуиция, не более.

[--mS--]

См. теорема 23 и следствие 4.

[Sinclair]

--mS-- в сообщении #588529 писал(а):

См. теорема 23 и следствие 4.

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