Поиск перекрытия нормальных распределений с использованием двойного интеграла (dblquad) в MATLAB. Странное поведение

Я вычисляю перекрытие двух нормальных двумерных распределений, используя следующую функцию

function [ oa ] = bivariate_overlap_integral(mu_x1,mu_y1,mu_x2,mu_y2)
    %calculating pdf. Using x as vector because of MATLAB requirements for integration
    bpdf_vec1=@(x,y,mu_x,mu_y)(exp(-((x-mu_x).^2)./2.-((y-mu_y)^2)/2)./(2*pi));

    %calcualting overlap of two distributions at the point x,y
    overlap_point = @(x,y) min(bpdf_vec1(x,y,mu_x1,mu_y1),bpdf_vec1(x,y,mu_x2,mu_y2));

    %calculating overall overlap area
    oa=dblquad(overlap_point,-100,100,-100,100);

Вы можете видеть, что это включает в себя получение двойного интеграла (x: -100 до 100, y: -100 до 100, в идеале -inf до inf, но достаточно на данный момент) от функции перекрытия_точки, которая составляет минимум 2 pdf-s, заданных функция bpdf_vec1 двух распределений в точке x,y.

Теперь PDF никогда не равен 0, поэтому я ожидаю, что чем больше площадь интервала, тем больше будет конечный результат, очевидно, с незначительной разницей после определенной точки. Однако оказывается, что иногда, когда я уменьшаю размер интервала, результат растет. Например:

>> mu_x1=0;mu_y1=0;mu_x2=5;mu_y2=0;
>> bpdf_vec1=@(x,y,mu_x,mu_y)(exp(-((x-mu_x).^2)./2.-((y-mu_y)^2)/2)./(2*pi));
>>  overlap_point = @(x,y) min(bpdf_vec1(x,y,mu_x1,mu_y1),bpdf_vec1(x,y,mu_x2,mu_y2));
>> dblquad(overlap_point,-10,10,-10,10)
ans =
    0.0124
>> dblquad(overlap_point,-100,100,-100,100)
ans =
    1.4976e-005  -----> strange, as theoretically cannot be smaller then the first answer
>> dblquad(overlap_point,-3,3,-3,3)
ans =
    0.0110  -----> makes sense that the result is less than the first answer as the           
interval is decreased

Здесь мы можем проверить, что перекрытия равны (близки) к 0 в граничных точках интервала.

>> overlap_point (100,100)
ans =
 0
>> overlap_point (-100,100)
ans =
 0
>> overlap_point (-100,-100)
ans =
 0
>> overlap_point (100,-100)
ans =
 0

Возможно, это связано с реализацией dblquad, или я где-то ошибаюсь? Я использую MATLAB R2011a.

Спасибо


matlab statistics integral normal-distribution
person shiftyscales    schedule 16.04.2013    source источник
comment
Я проверил ваш код и не вижу ошибок, поэтому я предполагаю, что это проблема точности с Matlab (точность каким-то образом связана с диапазоном ограничений). Я протестирую ваш код в октаве и посмотрю, есть ли у него такая же проблема. Кстати, вы можете попробовать это с необязательным параметром точности, например, например, dblquad (-100,100,-100,100,1e-8).   -  person Stuart    schedule 16.04.2013
comment
Интересно, что при запуске точно такого же кода в gnu-Octave проблема не возникала. С точностью по умолчанию практически не было разницы в результате интеграла по [-10,10,-10,10] по сравнению с интегралом по [-100,100,-100,100]. Оба результата были около 0,01241933 btw.   -  person Stuart    schedule 16.04.2013
comment
Можете ли вы протестировать, постепенно увеличивая пределы, и посмотреть, где он впервые начинает давать сбой. Например, [-30,30,-30,30], затем [-40,40,-40,40] и так далее, пока результаты не начнут расходиться. Затем в этот момент попробуйте добавить необязательный параметр точности 1e-7 или 1e-8 и т. д. и посмотрите, решит ли это проблему. КСТАТИ. Моя версия Октавы 3.6.1   -  person Stuart    schedule 16.04.2013
comment
Спасибо за комментарии. Результат начинает уменьшаться с [-10,10,-10,10] >> dblquad(overlap_point,-9,9,-9,9) ans = 0.012434302954144 >> dblquad(overlap_point,-10,10,-10,10) ans = 0.012426083045525. Указание (достаточно большого) параметра точности, похоже, имеет положительный эффект, но не решает проблему полностью, например: >> dblquad(-9,9,-9,9,1e-10) ans = 0.012419331486452, dblquad(-10,10,-10,10,1e-10) ans = 0.012419330851169, dblquad(-100,100,-100,100,1e-10) ans = 0.012415587111195 Результаты намного ближе друг к другу, но все же уменьшаются от [-10,10, -10,10].   -  person shiftyscales    schedule 16.04.2013
comment
Эти значения в октаве немного уменьшаются или увеличиваются?   -  person shiftyscales    schedule 16.04.2013
comment
Эти значения в Octave немного уменьшаются или увеличиваются? С точностью по умолчанию результат немного уменьшается при переходе от [-10,10,-10,10] к [-100,100,-100,100], поэтому да, это нелогично, но разница невелика (около -4e-12, то есть меньше целевой точности). Однако, когда я переоценил их, оба с точностью 1e-8, то результат [-100,100,-100,100] был на самом деле большим из двух (но только около 1e-12).   -  person Stuart    schedule 16.04.2013

Ответы (2)


arrow_upward
1
arrow_downward

ПОЗДРАВЛЯЕМ! Вы получаете награду за то, что стали 12-миллионным человеком, задавшим по сути этот вопрос. :) Я пытаюсь сказать, что это проблема, о которую все спотыкаются поначалу. Честно говоря, этот вопрос задают снова и снова, так что на самом деле этот вопрос следует пометить как дубликат.

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

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

Кстати, это НЕ то, что я ожидал увидеть в Octave и MATLAB. Это только проблема адаптивного интегратора и того, где он решит установить свои точки выборки.

person Community    schedule 16.04.2013
comment
так что, по сути, единственное решение имеет меньшие границы для интеграла? - person shiftyscales; 16.04.2013
comment
На самом деле я не жаловался, просто сказал, что вы не одиноки в этой проблеме. Я предполагаю, что мы все сделали в тот или иной момент. Тщательный, разумный выбор границ ДЕЙСТВИТЕЛЬНО поможет. Однако если вы станете слишком маленьким, вы потеряете часть массы. Некоторые квадратурные инструменты позволяют указать точки интереса, поэтому, если вы сможете указать пики нормалей как таковые, вы будете далеко впереди. Я не припоминаю, чтобы у dblquad была такая возможность. Я написал двухмерные инструменты адаптивной интеграции, и эту проблему никогда не бывает просто решить автоматически. - person ; 16.04.2013

arrow_upward
1
arrow_downward

Хорошо, вот мои результаты по октаве.

>fomat long
>z10 = dblquad(overlap_point,-10,10,-10,10)
z10 =  0.0124193303284589
>z100 = dblquad(overlap_point,-100,100,-100,100)
z100 =  0.0124193303245106
>z100 - z10
ans = -3.94835379669001e-012

>z10a = dblquad(overlap_point,-10,10,-10,10,1e-8)
z10a =  0.0124193306514996
>z100a = dblquad(overlap_point,-100,100,-100,100,1e-8)
z100a =  0.0124193306527155
>z100a-z10a
ans = 1.21588676627038e-012

КСТАТИ. Я замечал этот тип проблемы раньше с численными решениями. Иногда вы вносите изменение, которое, как вы ожидаете, повысит точность результата (в данном случае приближая пределы к идеальному случаю полной плоскости), но вместо этого вы получаете противоположный эффект, и результат становится менее точным. Что происходит в этом случае, так это то, что, выходя «шире», на -100..100, вы смещаете фокус с того места, где в вашей функции происходит действительно важное действие, которое близко к началу координат. В какой-то момент реализация dblquad, которую вы используете, должна начать увеличивать межсемпловое расстояние по мере увеличения пределов, а затем она начинает упускать некоторые важные вещи вблизи начала координат.

Может быть, кто-то с более поздней версией Matlab может проверить это и посмотреть, было ли оно улучшено.

person Stuart    schedule 16.04.2013