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 r).
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: np.sum(((v_model-v_obs)/errors)**2) where errors is an array of the error at each different r value. v_obs is the observed velocity at each r value (r is just an array of numbers). To minimize the function, I have to manipulate v_model which can be done via manipulating two "fixing parameters" p0 and r0 in the equation (integrand is shown below):
np.sqrt(4.302*10**(-6)*quad(integrand,0,r.all(),args=(p0,r0))[0]/r.all())
Before I get into the problem, I want to know if r.all() is appropriate since it wouldn't allow me to put r due to it being an array. An alternative I had was to make an array of v_model via:
#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))
Regardless, to minimize the function np.sum(((v_model-v_obs)/errors)**2)
I tried to do something like this:
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)
However, the values I get back aren't good fits at all (which is evident when I graph the observed data and my model)
In conclusion, I have two main question:
1.) Is my function taking the model minus the observed at each different r value, and if not, how do I fix this? (I think I messed up my v_model equation)
2.) Why is it returning wrong numbers for r0 and p0?
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)
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)
Let me know if I left out any data and thank you in advance!
If you change p0 in the function it brings the number down in the optimizer routine. That seems to work well, I guess the numbers are so large it causes some numerical errors in the solver. Just remember to multiply the p0 answer by 3.85e+09, just like is done in the chisqfunc routine.
from scipy.optimize import*
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate
from scipy.integrate import quad
#-----------------------------------------------------------------------------
am = 18 # number of measurements
#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
])
grav_const = 4.302e-06
#print "np.pi, grav_const = " + str(np.pi) + ", " + str(grav_const)
#-----------------------------------------------------------------------------
def func_to_integrate(rx, p0, r0):
rho = (p0 * r0**3) / ( (rx + r0) * (rx**2 + r0**2) )
function_result = rho * 4.0 * np.pi * rx**2
#print "rx, p0, r0 = " + str(rx) + ", " + str(p0) + ", " + str(r0)
return function_result
#-----------------------------------------------------------------------------
def simpsonsRule(func, a, b, n, p0, r0):
if n%2 == 1:
return "Not applicable"
else:
h = (b - a) / float(n)
s = func(a, p0, r0) + sum((4 if i%2 == 1 else 2) * func(a+i*h, p0, r0) for i in range(1,n)) + func(b, p0, r0)
return s*h/3.0
#-----------------------------------------------------------------------------
def chisqfunc(iter_vars):
global v_model
p0 = iter_vars[0] * 3.85e+09
r0 = iter_vars[1]
v_model = []
for index in range(0, am):
integral_result = simpsonsRule(func_to_integrate, 0.0, r[index], 200, p0, r0)
v_mod_result = np.sqrt(grav_const * integral_result / r[index])
v_model.append(v_mod_result)
chisq = 0.0
for index in range(0, am):
chisq += ((v_model[index] - v_obs[index])/errors[index])**2
#chisq += ((v_model[index] - v_obs[index]))**2
chisq = np.sqrt(chisq)
#chisq = np.sum(((v_model - v_obs)/errors)**2)
print "chisq = " + str(chisq)
#print "iterated p0, r0 = " + str(p0) + ", " + str(r0)
return chisq
#-----------------------------------------------------------------------------
initial_guess = np.array([1.0, 2.0])
#initial_guess = np.array([3.5*10**8.62, 1.5])
#resolution = minimize(chisqfunc, initial_guess, method='TNC') # , tol=1e-12, options={'accuracy':0})
#resolution = minimize(chisqfunc, initial_guess, method='TNC', options={'rescale':-1})
#resolution = minimize(chisqfunc, initial_guess, method='CG')
resolution = minimize(chisqfunc, initial_guess, method='Nelder-Mead')
# # options tried
# 'gtol':1e-6, 'xtol':-1, 'rescale':0, 'ftol':0, 'accuracy':0, 'offset':None, 'scale':None
print("This is the function",resolution)
#print "resolution.x"
#print resolution.x
print "iterated p0 = " + str(resolution.x[0])
print "iterated p0 * 3.85e+09 = " + str(resolution.x[0] * 3.85e+09)
print "iterated r0 = " + str(resolution.x[1])
print "p0 should be about = " + str(3.5e+08)
print "r0 should be about = " + str(1.5)
print "-----------------------------------------------"
print "iterated model errors"
for index in range(0,am):
print "v_obs, v_model, error = " + str(v_obs[index]) + ", " + str(v_model[index]) + ", " + str(v_obs[index] - v_model[index])
chisq = 0.0
for index in range(0, am):
chisq += ((v_model[index] - v_obs[index])/errors[index])**2
#chisq += ((v_model[index] - v_obs[index]))**2
chisq = np.sqrt(chisq)
print "iterated chisq = " + str(chisq)
print ""
#-----------------------------------------------------------------------------
print "model errors using p0 = 3.5*10**8.62, r0 = 1.8"
p0 = 3.5*10**8.62
r0 = 1.8
# chisq = 118 for p0 = 3.5*10**8.62, r0 = 1.8
# chisq = 564 for p0 = 3.5*10**8 , r0 = 1.5
chisq = 0.0
for index in range(0, am):
integral_result = simpsonsRule(func_to_integrate, 0.0, r[index], 200, p0, r0)
v_mod = np.sqrt(grav_const * integral_result / r[index])
print "v_obs, v_mod, error = " + str(v_obs[index]) + ", " + str(v_mod) + ", " + str(v_obs[index] - v_mod)
#chisq += ((v_mod - v_obs[index]))**2
chisq += ((v_mod - v_obs[index])/errors[index])**2
chisq = np.sqrt(chisq)
print "using p0 = 3.5*10**8.62, r0 = 1.8, the chisq = " + str(chisq)
#-----------------------------------------------------------------------------