Использование scipy.optimize.minimize с дискретным ядром Гаусса

Question

Я использую scipy.optimize.minimize, чтобы попытаться определить оптимальные параметры функции плотности вероятности (PDF). Мой PDF включает в себя дискретное ядро ​​Гаусса ( https://en.wikipedia.org/wiki/Gaussian_function и https://en.wikipedia.org/wiki/Scale_space_implementation).

Теоретически, я знаю среднее значение PDF (где PDF должен быть в центре). Так что, если бы я рассчитал ожидаемое значение моего PDF, я должен восстановить среднее значение, которое я уже знаю. В моем PDF-файле выбраны дискретные значения n (которые никогда не должны быть отрицательными и должны начинаться с 0, чтобы иметь какой-либо физический смысл), и я пытаюсь определить оптимальное значение t ("коэффициент масштабирования"), чтобы восстановить среднее значение из PDF (который опять же я знаю заранее).

Мой минимальный рабочий пример для определения оптимального "коэффициента масштабирования" t следующий:

#!/usr/bin/env python3                                                                                                                                                                                      

import numpy as np
from scipy.special import iv
from scipy.optimize import minimize

def discrete_gaussian_kernel(t, n):
    return np.exp(-t) * iv(n, t)

def expectation_value(t, average):

    # One constraint is that the starting value                                                                                                                                                             
    # of the range over which I sample the PDF                                                                                                                                                              
    # should be 0.                                                                                                                                                                                          

    # Method 1 - This seems to give good, consistent results                                                                                                                                                
    int_average = int(average)
    ceiling_average = int(np.ceil(average))
    N = range(int_average - ceiling_average + 1,
              int_average + ceiling_average + 2)

    # Method 2 - The multiplicative factor for 'end' is arbitrary.                                                                                                                                          
    #            I should in principle be able make end be as large as                                                                                                                                      
    #            I want since the PDF goes to zero for large values of n,                                                                                                                                   
    #            but this seems to impact the result and I do now know why.                                                                                                                                 
    #start = 0                                                                                                                                                                                              
    #end = 2 * int(average)                                                                                                                                                                                 
    #N = range(start, end)                                                                                                                                                                                  

    return np.sum([n * discrete_gaussian_kernel(t, n - average) for n in N])

def minimize_function(t, average):
    return average - expectation_value(t, average)

if __name__ == '__main__':

    average = 8.33342
    #average = 7.33342                                                                                                                                                                                      

    solution = minimize(fun = minimize_function,
                        x0 = 1,
                        args = average)

    print(solution)
    t = solution.x[0]
    print('          solution t =', t)
    print('       given average =', average)
    print('recalculated average =', expectation_value(t, average))

У меня две проблемы с моим минимальным рабочим примером:

1) Код работает нормально для некоторых значений, которые я выбираю для переменной "среднее". Одним из примеров этого является значение 8,33342. Однако код не работает для других значений, например, 7.33342. В этом случае я получаю

 RuntimeWarning: overflow encountered in exp

поэтому я подумал, что, возможно, scipy.optimize.minimize выбирает неверное значение для t (например, большое отрицательное число). Я уверен, что это проблема, так как я напечатал значение t в функции hopeation_value, и t становится все более отрицательным. Поэтому я хотел бы добавить границы возможных значений того, что может принимать "t" ("t" не должно быть отрицательным). Глядя на документацию scipy.optimize.minimize, есть аргумент ключевого слова bounds. Итак, я попробовал:

solution = minimize(fun = minimize_function,
                    x0 = 1,
                    args = average,
                    bounds = ((0, None)))

но я получаю ошибку:

 ValueError: length of x0 != length of bounds

Я искал эту ошибку в stackru, и есть некоторые другие потоки, но я не нашел ничего полезного. Как я могу успешно установить границу?

2) Мой другой вопрос связан с тем, что scipy.optimize.minimize чувствителен к диапазону, в котором я вычисляю значение ожидания. Для среднего значения

 average = 8.33342

и метод расчета диапазона как

# Method 1 - This seems to give good, consistent results                                                                                                                                                
int_average = int(average)
ceiling_average = int(np.ceil(average))
N = range(int_average - ceiling_average + 1,
          int_average + ceiling_average + 2)

"пересчитанное среднее" составляет 8,3329696426. Но для другого метода (который имеет очень похожий диапазон),

# Method 2 - The multiplicative factor for 'end' is arbitrary.                                                                                                                                          
#            I should in principle be able make end be as large as                                                                                                                                      
#            I want since the PDF goes to zero for large values of n,                                                                                                                                   
#            but this seems to impact the result and I do now know why.                                                                                                                                 
start = 0
end = 2 * int(average)
N = range(start, end)

"пересчитанное среднее" составляет 8,31991111857. Диапазоны в каждом случае одинаковы, поэтому я не знаю, почему происходит такое большое изменение, особенно с учетом того, что мое пересчитанное среднее должно быть как можно ближе к истинному среднему. И если бы я расширил диапазон до больших значений (что я считаю разумным, так как там PDF сводится к нулю),

start = 0
end = 4 * int(average)
N = range(start, end)

"пересчитанное среднее" составляет 9.12939372912, что еще хуже. Так есть ли последовательный метод для расчета диапазона, чтобы восстановленное среднее всегда было как можно ближе к истинному среднему? Коэффициент масштабирования может принимать любое значение, поэтому я думаю, что scipy.optimize.minimize должен быть в состоянии найти коэффициент масштабирования для точного возврата к истинному среднему значению.