численное взятие интеграла

[AndreyL]

Такая задача: есть две группы случайных величин с многомерным нормальным распределением с разными параметрами (каждая группа со своими параметрами). Вопрос о принадлежности какого-то измерения той или иной группе решается по величине плотности распределения. Тогда ошибка классификации (вероятность неправильной классификации) может быть посчитана как $\int \min (P_1(\bar{x}), P_2(\bar{x})) \, d \bar{x}$, где $P_1(\bar{x}) $ и $P_2(\bar{x}) $ - плотности вероятности первой и второй группы соответственно.
Вопрос такой: как посчитать эту вероятность при больших размерностях пространства? Вольфрамовская Математика ругается и говорит, что интеграл ноль.

[Vince Diesel]

Можно попробовать метод Монте-Карло:

NIntegrate[..., Method->"MonteCarlo"]

или

NIntegrate[..., Method->"AdaptiveMonteCarlo"]

[AndreyL]

Так то да, но чего-то не так получается:
генерируем выборку, так чтобы не все корреляции были нулевыми

ClearAll[r]
seq = Join[Sequence @@ Table[Transpose[RandomReal[MultinormalDistribution[{0, 0}, {{1, r}, {r, 1}}],100]], {r, -0.8, 0.8, 0.4}]] // Transpose;
{n, m} = Dimensions[seq]

оцениваем по ней распределение, вектор средних и ковариационную матрицу:

ClearAll[a, s]
A = Array[a, m];
S = Array[s, {m, m}];
dist = EstimatedDistribution[seq, MultinormalDistribution[A, S]];
av = Mean[dist];
cov = Covariance[dist];

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

fu = PDF[dist, A]
V = Table[{a[j], av[[j]] - 3 Sqrt[cov[[j, j]]], av[[j]] + 3 Sqrt[cov[[j, j]]]}, {j, m}]
NIntegrate[Evaluate[fu], Evaluate[Sequence @@ V], Method -> "MonteCarlo"]
NIntegrate[Evaluate[fu], Evaluate[Sequence @@ V], Method -> "AdaptiveMonteCarlo"]

Для MonteCarlo у меня стабильно получается цифра в районе 0.8-0.9, для AdaptiveMonteCarlo в районе 0.4, а должно быть примерно 1. Причем расширение пределов интегрирования только ухудшает ситуацию.

[Vince Diesel]

Может, точек маловато. Тогда надо увеличивать их число: MaxPoints->10^4 или сколько надо.

[AndreyL]

Vince Diesel в сообщении #575171 писал(а):

Может, точек маловато. Тогда надо увеличивать их число: MaxPoints->10^4 или сколько надо.

Надо, оказывается, много. Поставил 10^9 (больше не берет), два интеграла считаются ночь (потому и молчал долго), результаты такие:
0.899727
1.10749
1.0164
0.559048
0.966356
0.986756
по идее должны быть все единицы. Постоянно говорит, что "The integral failed to converge after 1000000100 integrand evaluations."
Метод "MonteCarlo", поскольку "AdaptiveMonteCarlo" совсем чего-то непонятное дает.

Есть ли какая-то возможность взятия таких интегралов на современной технике?

[Vince Diesel]

Если размерность большая, то не знаю, что, кроме Монте-Карло, работало бы. Разве что символьно, что тоже вряд ли.

Ну, получается, что так можно получить точность порядка нескольких процентов. Вопрос только в том, устраивает ли это. А если надо больше точек, то можно брать среднее арифметическое. Если есть несколько процессоров, можно параллельно запускать.

[AndreyL]

Vince Diesel в сообщении #576815 писал(а):

Если размерность большая, то не знаю, что, кроме Монте-Карло, работало бы. Разве что символьно, что тоже вряд ли.

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

Vince Diesel в сообщении #576815 писал(а):

Ну, получается, что так можно получить точность порядка нескольких процентов. Вопрос только в том, устраивает ли это. А если надо больше точек, то можно брать среднее арифметическое. Если есть несколько процессоров, можно параллельно запускать.

Сложность в том, что эти несколько процентов и определяют задачу - здесь считается вероятность ошибки классификации. Вопрос - а что значит "брать среднее арифметическое"? Просчитать интеграл несколько раз и взять среднее? И как в Математике 8 запустить NIntegrate параллельно?

[Vince Diesel]

AndreyL в сообщении #577527 писал(а):

Просчитать интеграл несколько раз и взять среднее?

Да.

AndreyL в сообщении #577527 писал(а):

И как в Математике 8 запустить NIntegrate параллельно?

results = ParallelTable[NIntegrate[...], {k, 1, 12}]
Mean[results]

Число раз лучше брать пропорционально количеству ядер, чтобы все были загружены.

[AndreyL]

Спасибо! Буду пробовать.