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

,
где:

- число измерений.

- параметер масштабирования (scale).

- радиус-вектор в

-мерном пространстве.

- скалярное произведение

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

зависимы, то мы не можем сгенерировать вектор покомпонентно.
Однако, можно заметить, что для всех векторов одинаковой длины (то есть, векторов, у которых совпадает

), функция дает одинаковую плотность вероятности. Таким образом, в полярных координатах, плотность вероятности зависит только от длины

, и никак не зависит от направления.
Это позволяет нам разбить генерацию вектора на 2 части:
1) Генерацию направления вектора, что тривиально сводится к генерации равномерно распределенных углов для полярных координат, и не представляет дальнейшего интереса.
2) Генерацию длины вектора.
Генерацию длины вектора я полагаю осуществить следующим образом:
2a) Получить фукнцию плотности вероятности для длины вектора.
2b) Получить функцию распределения вероятности, проинтегрировав плотность вероятности.
2c) Найти обратную функцию к функции распределения вероятности.
Тогда, генерировать длину вектора можно будет просто выбрав равномерно распределенную случайную величину из области значений функции распределения вероятности, и, получив длину с помощью обратной функции.
Теперь подробнее:
2a) Нужно перейти от плотности вероятности для вектора, к плотности вероятности для длины вектора.
Как мне подсказывает математическая интуиция, и соображения о пределе отношения объема полых гиперсфер с большим и малым радиусом, плотность вероятности длины вектора получается домножением исходной плотности вероятности в произвольной точке с данной длиной на площадь поверхности гиперсферы с радиусом, равным этой длине (иначе говоря - на интеграл по всей области с данной длиной).
Языком формул, получается следующий результат (поскольку, мне не нужен нормировочный коэффициент, от радиуса сферы остается только

):

К сожалению, мои математические познания не достаточны для того, чтобы твердо доказать этот способ, и я даже не знаю, как это можно проверить.
2b) Теперь, имея плотность вероятности, я могу найти функцию распределения вероятности длины вектора.
Для этого мне нужно взять интеграл

. Сам взять подобный интеграл я вряд ли смогу, я взял его с помощью 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}$ $\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}$](https://dxdy-01.korotkov.co.uk/f/4/2/3/423714ed4a11f5367bbe9e7aa385eb3682.png)
,
где

- гипергеометрическая функция. С помощью, опять же, вольфрама, я могу найти предел этой функции при

, стремящимся к бесконечности, для определеня границ генерации равномерно распределенного случайного значения.
2c) После этого, остается найти обратную функцию, и генератор готов.
Здесь вознакают серьезные трудности, потому что нужно найти обратную фукнкцию к произведению показательной и гипергеометрической функций. Как это сделать я не представляю. К счастью, этого можно избежать. Первый способ - построить обратную фукнцию с помощью кусочно-линейной аппроксимации. Второй - используя то, что функция распределения монотонно возрастает, каждый раз искать обратное значение с помощью двоичного поиска.
Отсюда возникает ряд вопросов.
Первый и самый главный - правильно ли я вообще все делаю?
То есть, нет ли у меня структурной ошибки в написании генератора?
Если вы разбираетесь в теме, и видите, что все правильно, пожалуйста, напишите это явно.
Второй - можно ли, и если да, то как, найти в общем виде обратную функцию на шаге 2c? Не обязательно точное значение, подойдет также ряд, или любой способ получить решение с задаваемой точностью. Очевидно, что это способ должен быть быстрее двоичного поиска, иначе в нем нет смысла.
Третий. Я еще сам в этом направлении не копал, но задам вопрос, раз уж пишу - какие вы могли бы порекомендовать библиотеки для вычисления гипергеометрической функции? (ответ имеет смысл только если нельзя найти обратную функцию в общем виде без гипергеометрической). Пишу на Cпп под линукс, главное - скорость.