title: Using scipy.optimize.minimize with a discrete Gaussian kernel
body: <p>I am using scipy.optimize.minimize to try to determine the optimal parameters of a probability density function (PDF). My PDF involves a discrete Gaussian kernel (<a href="https://en.wikipedia.org/wiki/Gaussian_function" rel="nofollow noreferrer">https://en.wikipedia.org/wiki/Gaussian_function</a> and <a href="https://en.wikipedia.org/wiki/Scale_space_implementation#The_discrete_Gaussian_kernel" rel="nofollow noreferrer">https://en.wikipedia.org/wiki/Scale_space_implementation#The_discrete_Gaussian_kernel</a>).</p>
  
  <p>In theory, I know the average value of the PDF (where the PDF should be centered on). So if I were to calculate the expectation value of my PDF, I should recover the mean value that I already know. My PDF is sampled at discrete values of n (which must never be negative and should start at 0 to make any physical sense), and I am trying to determine the optimal value of t (the "scaling factor") to recover the average value of the PDF (which again, I already know ahead of time).</p>
  
  <p>My minimal working example to determine the optimal "scaling factor" t is the following:</p>
  
  <pre><code>#!/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))
  </code></pre>
  
  <p>I have two problems with my minimal working example:</p>
  
  <p>1) The code works OK for some values of what I choose for the variable "average." One example of this is when the value is 8.33342. However, the code does not work for other values, for example 7.33342. In this case, I get</p>
  
  <pre><code> RuntimeWarning: overflow encountered in exp
  </code></pre>
  
  <p>so I was thinking that maybe scipy.optimize.minimize is choosing a bad value for t (like a large negative number). I am confident that this is the problem since I have printed out the value of t in the function expectation_value, and t becomes increasingly negative. So I would like to add bounds to the possible values of what "t" can take ("t" should not be negative). Looking at the documentation of scipy.optimize.minimize, there is a bounds keyword argument. So I tried:</p>
  
  <pre><code>solution = minimize(fun = minimize_function,
                      x0 = 1,
                      args = average,
                      bounds = ((0, None)))
  </code></pre>
  
  <p>but I get the error:</p>
  
  <pre><code> ValueError: length of x0 != length of bounds
  </code></pre>
  
  <p>I searched for this error on stackoverflow, and there are some other threads, but I did not find any helpful. How can I set a bound successfully?</p>
  
  <p>2) My other question has to do with scipy.optimize.minimize being sensitive to the range over which I calculate the expectation value. For an average value of </p>
  
  <pre><code> average = 8.33342
  </code></pre>
  
  <p>and the method of calculating the range as</p>
  
  <pre><code># 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)
  </code></pre>
  
  <p>the "recalculated average" is 8.3329696426. But for the other method (which has a very similar range),</p>
  
  <pre><code># 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)
  </code></pre>
  
  <p>the "recalculated average" is 8.31991111857. The ranges are similar in each case, so I don't know why there is such a large change, especially since I what my recalculated average to be as close as possible to the true average. And if I were to extend the range to larger values (which I think is reasonable since the PDF goes to zero there),</p>
  
  <pre><code>start = 0
  end = 4 * int(average)
  N = range(start, end)
  </code></pre>
  
  <p>the "recalculated average" is 9.12939372912, which is even worse. So is there a consistent method to calculate the range so that the reconstructed average is always as close as possible to the true average? The scaling factor can take on any value so I would think scipy.optimize.minimize should be able to find a scaling factor to get back the true average exactly.</p>