衰减曲线最适合 SciPy
Decay curve best fit SciPy
我在尝试找到最适合我的数据时遇到了问题。使用 scipy.optimize.curve_fit 创建最佳匹配。我的数据和代码是：
编辑 您可以从here 下载数据文件。
数据是，
a b b2
55478 1.07E+43 54395.93833
56333 1.63E+43 54380.01385
57540 2.57E+43 52393.31605
61866 7.32E+43 52212.22838 52212.22838
代码：
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from __future__ import division
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import fit
import glob
import os
from scipy.optimize import curve_fit
import matplotlib.patches as patches
pf = pd.read_csv('/home/imhotep/Desktop/lala.csv', sep=',', encoding ='utf-8')
a1= pf['a'].max()
b1 = pf['b2'].max()
npoc=100
x = np.linspace((b1), (pf['b'].max()),npoc)
yy = np.linspace((pf['a'].min()), (pf['a'].max()), npoc)
fig = plt.figure()
ax4 = fig.add_subplot(111)
def h(x,k):
return a1* (((x-(b1))/(k))**(-(5./3.)))
popt,pcov = curve_fit(h,x,yy)
print 'POPT,', popt,'PCOV',pcov
y_fi1 = h(x, *popt)
ax4.plot(x, y_fi1, label='fit', ls='-', color='blue')
ax4.plot(pf['b'], pf['a'], ls='None', color='blue', marker='o')
plt.show()
就像那样。当我 运行 我得到合适的代码时：
但是，大致应该是这样的：
谁能告诉我哪里出错了？我是曲线拟合的初学者。
您想为 a 和 b 描述的 4 个蓝点拟合模型吗？
那你应该往这个方向做：
popt,pcov = curve_fit(h,b,a)
编辑：
如问题和此答案的评论中所述，您应该仅在原始数据上使用拟合函数，然后使用 np.linspace 显示拟合的新创建数组。
这是我从你的代码中得到的：
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from __future__ import division
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
pf = pd.read_csv('lala.csv', sep=',', encoding ='utf-8')
a1 = pf['a'].max()
#b1 = pf['b2'].max()
x = pf["b"]
y = pf["a"]
def h(x,k,b1):
return a1*((x-b1)/k)**(-5/3)
popt,pcov = curve_fit(h,x,y)
print 'POPT,', popt,'PCOV',pcov
xfit = np.linspace(x.min(),x.max(),100)
y_fi1 = h(xfit, *popt)
fig = plt.figure()
ax4 = fig.add_subplot(111)
ax4.plot(xfit, y_fi1, label='fit', ls='-', color='blue')
ax4.plot(x, y, ls='None', color='blue', marker='o')
plt.show()
使用 curve_fit 仅查找参数 k 会导致错误，因此我将 b1 作为搜索参数。然后它确实找到了合适的，但仍然不完全令人满意。
输出：
POPT, [ 238.09666313 51973.04601693]
PCOV [[ 21500.32886377 -22370.88448044] [-22370.88448044 23850.34961769]]
您可以尝试通过首先以 old-fashioned 方式线性化方程来获得更好的 k 初始估计...
...然后使用您选择的软件计算简单的线性回归。在这里，我使用 statsmodels 得到 0.0007 for 1/k，这意味着与 curve_fit.[= 一起使用的初始估计约为 1400 15=]
import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.api as sm
import matplotlib.pyplot as plt
df = pd.read_csv('lala.csv')
time_min = min(df.time)
luminosity_max = max(df.luminosity)
df['Y'] = (df.luminosity/luminosity_max)**(-0.6)
results = smf.ols('Y ~ time', data=df).fit()
print (results.summary())
fig, ax = plt.subplots()
fig = sm.graphics.plot_fit(results, 1, ax=ax)
plt.show()
从这段代码生成的图表中的误差条可以明显看出（如果不是这样的话）k 的估计存在相当大的不确定性。
我仍然没有成功 curve_fit 工作。但是，您可能对要优化的功能的这种修改感兴趣。重命名 csv 中的列后（这样变量对我来说就不那么混乱了）我以这种方式重写了主要代码。我发现 h 的值在 1400 左右第一次给出了 nan。我决定简单地用最大亮度替换这些 nan。如果你 运行 这个我想你会发现 k=700 给出了最好的（粗略的）开始。
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
df = pd.read_csv('lala.csv')
print (df)
luminosities = df.luminosity
times = df.time
luminosity_max= max(luminosities)
time_min = min(times)
def h(time,k):
result = luminosity_max*((time-time_min)/k)**(-5./3.)
if np.isnan(result):
result = luminosity_max
return result
for time in range(52000,56000,1000):
print (time, h(time,2800), h(time,1400),h(time,700))
#~ popt,pcov = curve_fit(h,times,luminosities,700)
#~ print ('POPT,', popt,'PCOV',pcov)
我不认为你在 curve_fit 上做错了什么。我怀疑数学模型非常不适合这些数据。这就是为什么。我 运行 以下代码计算 least-square 各种值 k 的误差。
import pandas as pd
import numpy as np
df = pd.read_csv('lala.csv')
print (df)
luminosities = df.luminosity
times = df.time
luminosity_max= max(luminosities)
time_min = min(times)
def h(time,k):
result = luminosity_max*(max(time-time_min,5)/k)**(-5./3.)
if np.isnan(result):
result = luminosity_max
return result
def LS(k):
return sum([(luminosity-h(time,k))**2 for (luminosity,time) in zip(luminosities,times)])
for k in range(10,110,10):
print (k, LS(k))
结果：
dummy1 luminosity time dummy2
0 55478 1.066349 54395.938333 NaN
1 56333 1.630938 54380.013854 NaN
2 57540 2.569603 52393.316048 NaN
3 61866 7.324060 52212.228380 52212.22838
10 263.810260704
20 4431.42454446
30 18991.1298817
40 51557.4862318
50 110648.507655
60 205525.705606
70 346090.508811
80 542810.685895
90 806664.876933
100 1149099.00104
我注意到k的值越小，LSE越小。但我认为这会使拟合模型 'hug' 成为水平轴，如您所见。
我在尝试找到最适合我的数据时遇到了问题。使用 scipy.optimize.curve_fit 创建最佳匹配。我的数据和代码是：
编辑 您可以从here 下载数据文件。 数据是，
a b b2
55478 1.07E+43 54395.93833
56333 1.63E+43 54380.01385
57540 2.57E+43 52393.31605
61866 7.32E+43 52212.22838 52212.22838
代码：
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from __future__ import division
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import fit
import glob
import os
from scipy.optimize import curve_fit
import matplotlib.patches as patches
pf = pd.read_csv('/home/imhotep/Desktop/lala.csv', sep=',', encoding ='utf-8')
a1= pf['a'].max()
b1 = pf['b2'].max()
npoc=100
x = np.linspace((b1), (pf['b'].max()),npoc)
yy = np.linspace((pf['a'].min()), (pf['a'].max()), npoc)
fig = plt.figure()
ax4 = fig.add_subplot(111)
def h(x,k):
return a1* (((x-(b1))/(k))**(-(5./3.)))
popt,pcov = curve_fit(h,x,yy)
print 'POPT,', popt,'PCOV',pcov
y_fi1 = h(x, *popt)
ax4.plot(x, y_fi1, label='fit', ls='-', color='blue')
ax4.plot(pf['b'], pf['a'], ls='None', color='blue', marker='o')
plt.show()
就像那样。当我 运行 我得到合适的代码时：
但是，大致应该是这样的：
谁能告诉我哪里出错了？我是曲线拟合的初学者。
您想为 a 和 b 描述的 4 个蓝点拟合模型吗？
那你应该往这个方向做：
popt,pcov = curve_fit(h,b,a)
编辑：
如问题和此答案的评论中所述，您应该仅在原始数据上使用拟合函数，然后使用 np.linspace 显示拟合的新创建数组。
这是我从你的代码中得到的：
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from __future__ import division
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
pf = pd.read_csv('lala.csv', sep=',', encoding ='utf-8')
a1 = pf['a'].max()
#b1 = pf['b2'].max()
x = pf["b"]
y = pf["a"]
def h(x,k,b1):
return a1*((x-b1)/k)**(-5/3)
popt,pcov = curve_fit(h,x,y)
print 'POPT,', popt,'PCOV',pcov
xfit = np.linspace(x.min(),x.max(),100)
y_fi1 = h(xfit, *popt)
fig = plt.figure()
ax4 = fig.add_subplot(111)
ax4.plot(xfit, y_fi1, label='fit', ls='-', color='blue')
ax4.plot(x, y, ls='None', color='blue', marker='o')
plt.show()
使用 curve_fit 仅查找参数 k 会导致错误，因此我将 b1 作为搜索参数。然后它确实找到了合适的，但仍然不完全令人满意。
输出：
POPT, [ 238.09666313 51973.04601693]
PCOV [[ 21500.32886377 -22370.88448044] [-22370.88448044 23850.34961769]]
您可以尝试通过首先以 old-fashioned 方式线性化方程来获得更好的 k 初始估计...
...然后使用您选择的软件计算简单的线性回归。在这里，我使用 statsmodels 得到 0.0007 for 1/k，这意味着与 curve_fit.[= 一起使用的初始估计约为 1400 15=]
import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.api as sm
import matplotlib.pyplot as plt
df = pd.read_csv('lala.csv')
time_min = min(df.time)
luminosity_max = max(df.luminosity)
df['Y'] = (df.luminosity/luminosity_max)**(-0.6)
results = smf.ols('Y ~ time', data=df).fit()
print (results.summary())
fig, ax = plt.subplots()
fig = sm.graphics.plot_fit(results, 1, ax=ax)
plt.show()
从这段代码生成的图表中的误差条可以明显看出（如果不是这样的话）k 的估计存在相当大的不确定性。
我仍然没有成功 curve_fit 工作。但是，您可能对要优化的功能的这种修改感兴趣。重命名 csv 中的列后（这样变量对我来说就不那么混乱了）我以这种方式重写了主要代码。我发现 h 的值在 1400 左右第一次给出了 nan。我决定简单地用最大亮度替换这些 nan。如果你 运行 这个我想你会发现 k=700 给出了最好的（粗略的）开始。
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
df = pd.read_csv('lala.csv')
print (df)
luminosities = df.luminosity
times = df.time
luminosity_max= max(luminosities)
time_min = min(times)
def h(time,k):
result = luminosity_max*((time-time_min)/k)**(-5./3.)
if np.isnan(result):
result = luminosity_max
return result
for time in range(52000,56000,1000):
print (time, h(time,2800), h(time,1400),h(time,700))
#~ popt,pcov = curve_fit(h,times,luminosities,700)
#~ print ('POPT,', popt,'PCOV',pcov)
我不认为你在 curve_fit 上做错了什么。我怀疑数学模型非常不适合这些数据。这就是为什么。我 运行 以下代码计算 least-square 各种值 k 的误差。
import pandas as pd
import numpy as np
df = pd.read_csv('lala.csv')
print (df)
luminosities = df.luminosity
times = df.time
luminosity_max= max(luminosities)
time_min = min(times)
def h(time,k):
result = luminosity_max*(max(time-time_min,5)/k)**(-5./3.)
if np.isnan(result):
result = luminosity_max
return result
def LS(k):
return sum([(luminosity-h(time,k))**2 for (luminosity,time) in zip(luminosities,times)])
for k in range(10,110,10):
print (k, LS(k))
结果：
dummy1 luminosity time dummy2
0 55478 1.066349 54395.938333 NaN
1 56333 1.630938 54380.013854 NaN
2 57540 2.569603 52393.316048 NaN
3 61866 7.324060 52212.228380 52212.22838
10 263.810260704
20 4431.42454446
30 18991.1298817
40 51557.4862318
50 110648.507655
60 205525.705606
70 346090.508811
80 542810.685895
90 806664.876933
100 1149099.00104
我注意到k的值越小，LSE越小。但我认为这会使拟合模型 'hug' 成为水平轴，如您所见。