title: Scipy.minimize not correctly minimizing
body: <p>I am currently doing a project that maps the velocity of stars as you move out from the center of the galaxy(the distance to the center is denoted by <code>r</code>).</p>
  
  <p>Essentially, I am aiming to minimize the distance between my model function and my observed function. To do this, I have to minimize the function: <code>np.sum(((v_model-v_obs)/errors)**2)</code> where <code>errors</code> is an array of the error at each different <code>r</code> value. <code>v_obs</code> is the observed velocity at each <code>r</code> value  (<code>r</code> is just an array of numbers). To minimize the function, I have to manipulate <code>v_model</code> which can be done via manipulating two "fixing parameters" <code>p0</code> and <code>r0</code> in the equation (integrand is shown below):</p>
  
  <pre><code>np.sqrt(4.302*10**(-6)*quad(integrand,0,r.all(),args=(p0,r0))[0]/r.all())
  </code></pre>
  
  <p>Before I get into the problem, I want to know if <code>r.all()</code> is appropriate since it wouldn't allow me to put <code>r</code> due to it being an array. An alternative I had was to make an array of <code>v_model</code> via:</p>
  
  <pre><code>#am is the amount of elements in the array r 
  #r, v_model,v_obs, and errors all have the same size
  def integrand(r,p0,r0):
      return (p0 * r0**3)/((r+r0)*((r**2)+(r0**2)))*4*3.1415926535*r**2
  integrals = []
  for i in r:
       integrals.append(quad(integrand, 0 ,i,args=(p0,r0)))
  
  v_model = []
  
  
  for x in range (0,am):
      k = integrals[x][0]
      i = r[x]
      v_model.append(np.sqrt((4.302*10**(-6)*k)/i))
  </code></pre>
  
  <p>Regardless, to minimize the function <code>np.sum(((v_model-v_obs)/errors)**2)</code>
  I tried to do something like this:</p>
  
  <pre><code>def chisqfunc(parameters):
      p0 = parameters[0]
      r0 = parameters[1]
      v_model = []
      for x in range(0,am):
          v_model.append(np.sqrt(4.302*10.0**(-6)*quad(integrand, 0, r[x], args=(p0,r0))[0]/r[x]))
      chisq = np.sum(((v_model-v_obs)/errors)**2)
      return chisq
  x0 = np.array([10**6,24])
  resolution = minimize(chisqfunc,x0)
  </code></pre>
  
  <p>However, the values I get back aren't good fits at all (which is evident when I graph the observed data and my model)</p>
  
  <p>In conclusion, I have two main question:</p>
  
  <p>1.) Is my function taking the model minus the observed at each different <code>r</code> value, and if not, how do I fix this? (I think I messed up my <code>v_model</code> equation)</p>
  
  <p>2.) Why is it returning wrong numbers for <code>r0</code> and <code>p0</code>?</p>
  
  <p>Here is my full code (By the way, to know if the minimization is working properly: r0 should be around 1.5 and p0 should be around: 3.5*10**8)</p>
  
  <pre><code>from scipy.optimize import*
  import numpy as np
  import matplotlib.pyplot as plt
  import scipy.integrate as integrate
  from scipy.integrate import quad
  
  #number of measurements
  am = 18
  
  
  r0 = 1.8
  p0 = 3.5*10**8.62
  
  #Observed Data
  v_obs = np.array([
  234.00,
  192.00,
  212.00,
  223.00,
  222.00,
  224.00,
  224.00,
  226.00,
  226.00,
  227.00,
  227.00,
  224.00,
  220.00,
  218.00,
  217.00,
  216.00,
  210.00,
  208.00
  
  
  
  ]
  )
  
  r = np.array([0.92, 
  2.32,
  3.24,   
  4.17,   
  5.10,   
  6.03,   
  6.96,   
  7.89,   
  8.80,   
  9.73,   
  10.64,
  11.58,
  12.52,
  13.46,
  14.40,
  15.33,
  16.27,
  17.11
  
  
  ]
  )
  
  errors = np.array([
  3.62,
  4.31,
  3.11,
  5.5,
  3.9,
  3.5,
  2.7,
  2.5,
  2.3,
  2.1,
  2.3,
  2.6,
  3.1,
  3.2,
  3.2,
  3.1,
  2.9,
  2.8
      ])
  
  
  #integral 
  def integrand(r,p0,r0):
      return (p0 * r0**3)/((r+r0)*((r**2)+(r0**2)))*4*3.1415926535*r**2
  integrals = []
  for i in r:
       integrals.append(quad(integrand, 0 ,i,args=(p0,r0)))
  
  v_model = []
  
  
  for x in range (0,am):
      k = integrals[x][0]
      i = r[x]
      v_model.append(np.sqrt((4.302*10**(-6)*k)/i))
  
  
  def chisqfunc(parameters):
      p0 = parameters[0]
      r0 = parameters[1]
      v_model = np.sqrt(4.302*10**(-6)*quad(integrand,0,r.all(),args=(p0,r0))[0]/r.all())
      chisq = np.sum(((v_model-v_obs)/errors)**2)
      print(v_model)
      return chisq
  x0 = np.array([10**6,24])
  resolution = minimize(chisqfunc,x0)
  print("This is the function",resolution)
  </code></pre>
  
  <p>Let me know if I left out any data and thank you in advance!</p>