1

当我试图找到最适合我的数据时，我遇到了问题。使用scipy.optimize.curve_fit创建最佳拟合。我的数据和代码是：

编辑您可以从这里下载数据文件。数据是， 

         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()

像那样。当我运行代码时，我变得很合适： 合身

但是，大致应该是这样的：

装的

谁能告诉我哪里出错了？我是曲线拟合的初学者。

4

3 回答 3

2

您想让模型适合 a 和 b 描述的 4 个蓝点吗？

那么你应该朝着这个方向努力：

popt,pcov = curve_fit(h,b,a)

编辑：

正如问题的评论和这个答案中提到的，您应该只对原始数据使用 fit 函数，然后使用新创建的数组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_fitto find only 参数会k导致错误，因此我将 b1 作为搜索参数。然后它确实找到了合适的位置，但仍然不能完全令人满意。输出：

POPT, [   238.09666313  51973.04601693] 
PCOV [[ 21500.32886377 -22370.88448044] [-22370.88448044  23850.34961769]]
于 2017-01-16T14:39:16.753 回答
0

k您可以尝试通过首先以老式方式线性化您的方程来获得更好的初始估计......

线性化形式

...然后使用您选择的软件计算简单的线性回归。在这里，我使用statsmodels得到 1/k 的 0.0007，这意味着与curve_fit一起使用的初始估计值约为 1400 。

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)
于 2017-01-17T00:09:04.047 回答
0

我不认为你在curve_fit上做错了什么。我怀疑数学模型非常不适合这些数据。这就是为什么。我运行以下代码来计算各种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 就越小。但我认为这会使拟合模型“拥抱”水平轴，如您所见。

于 2017-01-17T17:03:22.497 回答