scipy convolve зависит от x

Я пытаюсь свернуть логнормальный PDF и гауссовский PDF. Поэтому я определил функции следующим образом:

def PDF_log(x,sig,mu): # log normal PDF
   mu = np.log(mu)
   return( (1/x)*(1/(sig*np.sqrt(2*np.pi))) * np.exp(-(np.log(x)-mu)**2/(2*sig**2)) )   

def gauss(x,sig,mu): # a noraml PDF
   return( 1/(np.sqrt(2*np.pi*sig**2)) * np.exp(-(x-mu)**2/(2.*sig**2)) )

def gauss_log(x, sig, mu, sig0, mu0):
    a = signal.convolve(PDF_log(x,sig,mu),gauss(x,sig0,mu0),mode='same')/np.sum(gauss(x,sig0,mu0)) 


def test():
    mu = 0.6
    sig = 0.2
    sig0 = 0.05
    mu0 = mu
    x = np.linspace( 0.5, 0.6, 10000 )
    plt.plot( x, gauss_log(x, sig, mu, sig0, mu0), '--', label='gauss_X_log', zorder=10 )
    plt.plot( x, gauss(x,sig0,mu0), label='gauss' )
    plt.plot( x, PDF_log(x,sig,mu), label='log' )
    plt.legend()
    plt.show()

Это дает мне следующий результат:

первый домен массива

Красная линия — логнормальная PDF, зеленая — гауссова. «Свертка» — синяя пунктирная линия. Когда я меняю домен x x = np.linspace( 0.5, 0.8, 10000 ), я получаю совсем другие результаты:

для другого домена

Здесь явно что-то не так. Результат моего интеграла свертки F(x) = int (g(t)*f(x-t))dt не должен зависеть от диапазона "x". Затем я сделал домен большим, то есть x = np.linspace( 0.00001, 100, 10000 ), что дает мне такую ​​ерунду:

домен большого массива

Либо в моем сценарии есть простая ошибка, либо я неправильно понимаю дискретную свертку.


python scipy convolution
person Sebastiano1991    schedule 06.11.2017    source источник
comment
Вы можете сравнить свою логарифмическую норму с scipy.stats.lognorm   -  person percusse    schedule 06.11.2017
comment
Я думаю, что мой логарифмический PDF-файл в порядке, по крайней мере, я могу сопоставить их, когда использую scipy.stats.lognorm.pdf. Я просто скопировал формулу из Википедии. Я подозреваю, что моя проблема заключается в моем понимании того, что делает scipy.signal.convolve (синяя пунктирная линия на графиках).   -  person Sebastiano1991    schedule 06.11.2017
comment
@Sebastiano1991 Sebastiano1991, чтобы восстановить ожидаемые результаты, вы должны свернуть сигналы по всему их домену (от -бесконечности до бесконечности). Для нормальных распределений вы можете ограничить интеграцию средним значением - среднее значение в несколько раз больше стандартного отклонения + несколько раз в стандартное отклонение с разумной точностью, но здесь область значений слишком мала.   -  person Pierre de Buyl    schedule 06.11.2017
comment
Это то, чего я ожидал — поэтому я увеличил «область интеграции» на своем третьем графике до [~ 0,100]. Это все еще слишком мало для «меньшей, чем средняя сторона», то есть не будет ли это ошибкой для всех малых средних? Есть ли потенциальное обходное решение?   -  person Sebastiano1991    schedule 06.11.2017

Ответы (1)


arrow_upward
1
arrow_downward

Я понял, где у меня была ошибка:

вместо правильной функции ядра gauss(x0,sig0,mu); Я использовал тот же x для Гаусса.

Это делает то, что я считаю правильным (с PDF-файлами сверху):

def gauss_log(x, x0, sig, mu, sig0, mu0):
    a = signal.convolve(gauss(x0,sig0,mu0),PDF_log(x,sig,mu),mode='same')/np.sum(gauss(x0,sig0,mu0))
    return( a )

def test(lightcurve,noisecurve):
    mu = 0.1 #can now put arbitrary small values of mu
    sig = 0.1
    sig0 = 0.05
    mu0 = 0
    x = np.linspace( 0.00001, 5, 1000 )
    x0 = np.linspace(-5,5,1000) #note that arrays need to be equal length!
    g_log = gauss_log(x, x0, sig, mu, sig0, mu0)
    plt.plot( x, g_log, '--', label='gauss_X_log', zorder=10 )
    plt.plot( x0, gauss(x0,sig0,mu0), label='gauss' )
    plt.plot( x, PDF_log(x,sig,mu), label='log' )
    plt.legend()
    plt.show()

    ###testing normalization!
    print(np.trapz(gauss(x0,sig0,mu0),x0))
    print(np.trapz(g_log,x))
    print(np.trapz(PDF_log(x,sig,mu),x))

1.0

0.999938903253

1.0

введите здесь описание изображения

person Sebastiano1991    schedule 07.11.2017
comment
На примечании к сайту: scipy.signal.convolve — это дискретная свертка двух массивов. Для большого количества точек массива это работает отлично. В моем реальном случае, когда у меня есть только 11 элементов массива, это вызывает проблемы. Чтобы обойти это, я решил вычислить интеграл свертки численно с помощью scipy.integrate.quad, что значительно улучшает результаты за счет времени вычислений. В качестве альтернативы можно выбрать значения PDF, которые являются наиболее ограничивающими, например. вокруг медианы. - person Sebastiano1991; 16.11.2017