Матожидание модифицированной многомерной гауссианы.

[Sinclair]

Добрый день.

У меня есть следующая задача, которую я пытаюсь решить уже несколько дней:
Существует пространство $\mathbb R^k$.
На этом пространстве имеется функция плотности вероятности - $k$-мерное нормальное распределение $N_{initial} \equiv N_k(\mu, \Sigma)$, где $\mu$ - матожидание, а $\Sigma$ - матрица ковариации распределения.
Все компоненты вектора $\mu$ лежат в пределах от 0 до 1.
Число измерений $k$ может равняться десяткам или нескольким сотням.
Данное распределение модифицированно следующим образом:
1) Распределение домножено на каждый компонент вектора $X$ в некоторой натуральной степени, и разницу единицы и этого компонента, тоже в какой-то натуральной степени.
2) В любых точках, хотя бы одна координата которых меньше 0 или больше 1 значение функции равно нулю.
Более наглядно напишу формулой:
$X \in \mathbb R^k$ - вектор, для которого вычисляется значение функции плотности вероятности.
$N_{mod}$ - модифицированное распределение
$x_1, x_2, ... x_k \in \mathbb R$ - компоненты вектора $x$.
$n_1, n_2, ... n_k \in \mathbb N$ - степени компонентов $x$.
$m_1, m_2, ... m_k \in \mathbb N$ - степени разниц 1 и компонентов $x$.
Степени могут достигать десятков, и даже сотен.
1) $N_{mod}(X) = x_1^{n_1} \cdot (1 - x_1)^{m_1} \cdot x_2^{n_2} \cdot (1 - x_2)^{m_2} \cdot ... \cdot x_k^{n_k} \cdot (1 - x_k)^{m_k} \cdot N_{initial}(X)$
2) $N_{mod}(X) = 0$, если, для какого либо $i$, $x_i < 0$, или $x_i > 1$
Мне нужно найти математическое ожидание получившейся функции распределения. То есть, решить задачу - чему будет равно матожидание вектора $X$, если он распределен согласно фукнции плотности вероятности $N_{mod}$.

Аналитический подход - взять интергал $\int X \cdot N_{mod}(X) dX$. Учитывая, что функция весьма сложная, развести переменные по сомножителям нельзя, и нет интегралов в элементарных функциях - это вряд ли можно считать тривиальной задачей.
Я пытался развести переменные по разным сомножителям, чтобы взять интеграл, но если я ввожу замену переменных, которая делает матрицу ковариации диагональной, то множители полиномной части функции становятся зависимыми от нескольких переменных. Для того, чтобы перемножить все в один полином, и поодиночке взять интегралы вида $\int x^n \cdot e^{-x^2} dx$ (в данном случае, $x$ - не вектор, а обычная вещественная переменная) - нужно астрономически много времени, так как степени могут быть достаточно большие - десятки, или даже сотни. Если у кого-то есть идеи, как взять такой интеграл - я буду очень благодарен.
Но в целом, мне хватит и численного решения этой задачи. Однако, проблемы возникают и тут. Учитывая, что число измерений $k$ может быть достаточно большим, методы вычисления интегралов, вроде метода прямоугольника и его вариаций не подходят.
Остается Монте-Карло. Однако, специфика фукнции $N_{mod}$ такова, что у нее часто будет один узкий пик, а остальная часть функции будет иметь на десятки порядков меньшее значение. Для того, чтобы в этот пик попало много значений (чтобы точность была удовлетворительной (мне нужна точность порядка 0.1%)), мне придется моделировать гигантское число значений. Это весьма нежелательно вычислительно.
Возможно, я как-то смогу ускорить выполнение Монте-Карло, если буду знать, где находится пик. Возможно, следует в его окрестность провести более тщательную симуляцию. Пока я остановился на этом, но вариант далек от идеального. Так что, был бы рад услышать советы и по более эффективному численному решению этой задачи.

[_hum_]

А попробовать использовать представление

$$I_i = \int x_i N_{mod}(X) dX = \int x_i h(X)N_{initial}(X)dX = \int x_i N_{initial}(X) dH(x) = \mathbf{E}\big(\xi_i N_{initial}(\overset{\rightarrow}{\xi})\big),$$

где $\overset{\rightarrow}{\xi}$ - случайный вектор, имеющий распределение $dH(X) = h(X)dX $ (многомерное Бета распределение вроде).
Тогда в методе Монте-Карло нужно будет генерировать реализации случайного вектора $\overset{\rightarrow}{\xi}$. А поскольку он распределен на $[0,1]^k$, то, возможно, "холостых" значений будет меньше.

[Sinclair]

_hum_ в сообщении #558505 писал(а):

А попробовать использовать представление

Что такое представление? Это какая-то разновидность замены? Я что-то не помню такого из университетской программы. Впрочем, времени прошло много.
Я правильно понимаю, вы рекомендуете мне генерировать некоторый вектор, в соответствии с бета-распределением (спасибо, что сказали, что это такое), а потом использовать результаты такой генерации для нахождения матожидания от исходного распределения так, как если бы я искал его методом монте-карло по случайным векторам?
Спасибо, я подумаю над этим (сейчас уже голова не варит).
Сразу вопрос, как генерировать вектор по этому распределению. Впрочем, я поищу, может не сложно.
А вообще, если все так, как я понимаю, то действительно, это похоже на хороший вариант для монте-карло. В смысле, для "ловли" пиков. Пиков там, впрочем, два будет. Еще один на гауссиане. Впрочем, его, видимо, "ловить" можно так же.

[_hum_]

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

Что такое представление? Это какая-то разновидность замены? Я что-то не помню такого из университетской программы.

Ну просто возьмите с обратного конца цепочки равенств распишите, убедитесь.

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

Я правильно понимаю, вы рекомендуете мне генерировать некоторый вектор, в соответствии с бета-распределением (спасибо, что сказали, что это такое), а потом использовать результаты такой генерации для нахождения матожидания от исходного распределения так, как если бы я искал его методом монте-карло по случайным векторам?

Да.

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

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

Да, алгоритмы генерации есть в сети. Вот только первый попавшийся мне как-то не очень порадовал по cкорости.

P.S. Про нормирующую константу в вашем распределении не забудьте.

[Sinclair]

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

_hum_ в сообщении #558675 писал(а):

Да, алгоритмы генерации есть в сети. Вот только первый попавшийся мне как-то не очень порадовал по cкорости.

Ну, с бета-распределением все просто. Его существует достаточно ограниченное число разновидностей. Если, например, $n$ и $m$ ограничены сотней, то можно просто заранее нагенирировать длинный массив случайных чисел для каждого распределения (всего будет 10к массивов). А потом просто читать нужный. Дешево и сердито. А так как в бета-распределении переменные разведены по сомножителям - можно генерировать раздельно по каждому измерению.
Плохо то, что пика будет два, причем, второй будет тоже достаточно узким, из-за проклятия размерности. При небольшом отклонении от матожидания по многим измерениям значение будет падать очень быстро. В принципе получается та же проблема. А многомерное кореллированное нормальное распределение генерировать сложнее. На данный момент, самая большая моя проблема в этом. Впрочем, думаю, что и для нормального кореллированного многомерного распределения генераторы в сети есть. Я просто не искал, потому что занят в другой части кода, пока.
Ну еще раз спасибо, что помогли. Так и буду решать, видимо.