2D gaussian distribution не суммируется с одним?

Question

Я построил обернутое двумерное гауссовское распределение в Python, используя приведенное здесь уравнение: http://www.aos.wisc.edu/~dvimont/aos575/Handouts/bivariate_notes.pdf Однако я не понимаю, почему мое распределение не суммируется до 1, несмотря на включение константы нормализации.

Для U х U решетки,

import numpy as np 
from math import * 

U = 60 
m = np.arange(U) 
i = m.reshape(U,1) 
j = m.reshape(1,U) 

sigma = 0.1 
ii = np.minimum(i, U-i) 
jj = np.minimum(j, U-j) 
norm_constant = 1/(2*pi*sigma**2) 
xmu = (ii-0)/sigma; ymu = (jj-0)/sigma 
rhs = np.exp(-.5 * (xmu**2 + ymu**2)) 
ker = norm_constant * rhs 

>> ker.sum() # area of each grid is 1 
15.915494309189533 

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

UPDATE:

Благодаря содержательные предложения других, я переписал код для применения L1 нормализации ядра. Тем не менее, представляется, что в контексте 2D свертка с помощью FFT, сохраняя диапазон как [0, U] может еще вернуться убедительный результат:

U = 100 
Ukern = np.copy(U) 
#Ukern = 15 

m = np.arange(U) 
i = m.reshape(U,1) 
j = m.reshape(1,U) 

sigma = 2. 
ii = np.minimum(i, Ukern-i) 
jj = np.minimum(j, Ukern-j) 
xmu = (ii-0)/sigma; ymu = (jj-0)/sigma 
ker = np.exp(-.5 * (xmu**2 + ymu**2)) 
ker /= np.abs(ker).sum() 

''' Point Density ''' 
ido = np.random.randint(U, size=(10,2)).astype(np.int) 
og = np.zeros((U,U)) 
np.add.at(og, (ido[:,0], ido[:,1]), 1) 

''' Convolution via FFT and inverse-FFT ''' 
v1 = np.fft.fft2(ker) 
v2 = np.fft.fft2(og) 
v0 = np.fft.ifft2(v2*v1) 
dd = np.abs(v0) 

plt.plot(ido[:,1], ido[:,0], 'ko', alpha=.3) 
plt.imshow(dd, origin='origin') 
plt.show() 

С другой стороны, проклейки ядро используя закомментированные линии дает неверную сюжет:

Answer 1

ПРИМЕЧАНИЕ: Как указано в комментариях сильфонных, это решение действует только тогда, когда вы пытаетесь построить гауссово ядро ​​свертки (или гауссов фильтр) для обработки изображений. Это не правильно нормализованная функция плотности гауссова, но это форма, которая используется для удаления гауссовского шума из изображений.

---

Вы пропускаете нормализацию L1:

ker /= np.abs(ker).sum() 

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

Фактически, константа гауссовой норнализации, которую вы можете использовать, может быть использована только для использования нормы L1 выше. Если я не worng, вы пытаетесь создать гауссовую свертку, а th выше - обычная нормализация tecnique, применяемая к ней.

Ваша вторая ошибка, как заявил @Praveen, заключается в том, что вам нужно пробовать сетку от [-U//2, U//2]. Вы можете сделать это:

i, j = np.mgrid[-U//2:U//2+1, -U//2:U//2+1] 

Наконец, если то, что вы пытаетесь сделать, это построить гауссовский фильтр, размер ядра, как правило, оценивается от сигмы (чтобы избежать нулей далеко от центра), как U//2 <= t * sigma , где t является параметром усечения, обычно установленным t=3 или t=4.

Answer 2

В настоящее время (значительно увеличенный в) контурном график вашего ker выглядит следующим образом:

Как вы можете видеть, это не выглядит как гауссово ядро. Большая часть вашей функции отмирает от 0 до 1.Глядя на самом ядре показывает, что все значения действительно отмирают действительно очень быстро:

>>> ker[0:5, 0:5] 
array([[ 1.592e+001, 3.070e-021, 2.203e-086, 5.879e-195, 0.000e+000], 
     [ 3.070e-021, 5.921e-043, 4.248e-108, 1.134e-216, 0.000e+000], 
     [ 2.203e-086, 4.248e-108, 3.048e-173, 8.136e-282, 0.000e+000], 
     [ 5.879e-195, 1.134e-216, 8.136e-282, 0.000e+000, 0.000e+000], 
     [ 0.000e+000, 0.000e+000, 0.000e+000, 0.000e+000, 0.000e+000]]) 

Значение суммы 15.915, что вы получаете в основном только кекли [0, 0]. Что все это говорит вам, так это то, что вы неправильно строите свою сетку.

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

Итак, во-первых, если вы хотите полную плотность вокруг mu=0, вы должны принять i и j от -U // 2 до U // 2. Но для решения вашей проблемы с разрешением я рекомендую принимать U количество точек между -0,5 и 0,5.

import numpy as np 
import matplotlib.pyplot as plt 

U = 60 
m = np.linspace(-0.5, 0.5, U) # 60 points between -1 and 1 
delta = m[1] - m[0]    # delta^2 is the area of each grid cell 
(x, y) = np.meshgrid(m, m)  # Create the mesh 

sigma = 0.1 
norm_constant = 1/(2 * np.pi * sigma**2) 

rhs = np.exp(-.5 * (x**2 + y**2)/sigma**2) 
ker = norm_constant * rhs 
print(ker.sum() * delta**2) 

plt.contour(x, y, ker) 
plt.axis('equal') 
plt.show() 

Это дает сумму, близкую к 1,0, а ядро ​​с центром в mu=0, как и ожидалось.

Зная, какой диапазон выбора (от -0,5 до 0,5) в этом случае зависит от вашей функции. Например, если вы сейчас возьмете sigma = 2, вы обнаружите, что ваша сумма не сработает, потому что теперь вы отбираете слишком мелко. Настройка вашего диапазона будет функцией ваших параметров - что-то вроде -5 * sigma - 5 * sigma - может быть лучшим вариантом.