【发布时间】：2017-12-05 08:52:38
【问题描述】：
我目前正在做一个项目，绘制从星系中心向外移动时恒星的速度（到中心的距离用r 表示）。
本质上，我的目标是最小化我的模型函数和我观察到的函数之间的距离。为此，我必须最小化函数：np.sum(((v_model-v_obs)/errors)**2) 其中errors 是每个不同r 值处的错误数组。 v_obs 是在每个 r 值处观察到的速度（r 只是一个数字数组）。为了最小化函数，我必须操纵v_model，这可以通过操纵方程式中的两个“固定参数”p0 和r0 来完成（积分如下所示）：
np.sqrt(4.302*10**(-6)*quad(integrand,0,r.all(),args=(p0,r0))[0]/r.all())
在我解决这个问题之前，我想知道r.all() 是否合适，因为它不允许我放置r，因为它是一个数组。我的另一种选择是通过以下方式制作v_model 数组：
#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))
不管怎样，最小化函数np.sum(((v_model-v_obs)/errors)**2)
我试图做这样的事情：
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)
但是，我得到的值根本不合适（当我绘制观察到的数据和我的模型时，这一点很明显）
最后，我有两个主要问题：
1.) 我的函数是否将模型减去在每个不同的r 值处观察到的值，如果不是，我该如何解决这个问题？ （我想我搞砸了我的v_model 方程）
2.) 为什么返回 r0 和 p0 的错误数字？
这是我的完整代码（顺便说一下，要知道最小化是否正常工作：r0 应该在 1.5 左右，p0 应该在：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)
如果我遗漏了任何数据，请告诉我，并提前感谢您！
【问题讨论】：
-
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嗯...使用您的代码，您的 最佳 x 计算结果为 ~19k。有一个解决方案~134.79可能还有更好的解决方案。你自己的解释也是如此。
标签： python scipy minimization