使用 scipy 在符号图中进行曲线拟合或插值答案

【问题标题】：Curve fit or interpolation in a semilogy plot using scipy使用 scipy 在符号图中进行曲线拟合或插值
【发布时间】：2018-07-07 01:26:25
【问题描述】：

我的数据点很少，我想创建一条线以在以符号学比例绘制时最适合数据点。我已经尝试过 scipy 的曲线拟合和三次插值，但与数据趋势相比，它们对我来说似乎都不是很合理。

请您检查是否有更有效的方法来创建适合数据的直线。可能外推可以做，但我没有找到关于python外推的好文档。

非常感谢您的帮助

import sys
import os
import numpy
import matplotlib.pyplot as plt
from pylab import *
from scipy.optimize import curve_fit
import scipy.optimize as optimization
from scipy.interpolate import interp1d
from scipy import interpolate

Mass500 = numpy.array([ 13.938 , 13.816,  13.661,  13.683,  13.621,  13.547,  13.477,  13.492,  13.237,
  13.232,  13.07,   13.048,  12.945,  12.861,  12.827,  12.577,  12.518])

y500 = numpy.array([  7.65103978e-06,   4.79865790e-06,   2.08218909e-05,   4.98385924e-06,
   5.63462673e-06,   2.90785458e-06,   2.21166794e-05,   1.34501705e-06,
   6.26021870e-07,   6.62368879e-07,   6.46735547e-07,   3.68589447e-07,
   3.86209019e-07,   5.61293275e-07,   2.41428755e-07,   9.62491134e-08,
   2.36892162e-07])



plt.semilogy(Mass500, y500, 'o')


# interpolation
f2 = interp1d(Mass500, y500, kind='cubic')
plt.semilogy(Mass500, f2(Mass500), '--')


# curve-fit
def line(x, a, b):
    return 10**(a*x+b)

#Initial guess.
x0     = numpy.array([1.e-6, 1.e-6])

print optimization.curve_fit(line, Mass500, y500, x0)
popt, pcov = curve_fit(line, Mass500, y500)
print popt
plt.semilogy(Mass500, line(Mass500, popt[0], popt[1]), 'r-')



plt.legend(['data', 'cubic', 'curve-fit'], loc='best')

show()

【问题讨论】：

标签： python scipy interpolation curve-fitting extrapolation

【解决方案1】：

numpy 和 scipy 中有许多回归函数可用。 scipy.stats.lingress 是比较简单的函数之一，它返回常见的线性回归参数。

这里有两个拟合半对数数据的选项：

绘制转换后的数据
重新缩放轴并转换输入/输出函数值

给定

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

%matplotlib inline


# Data
mass500 = np.array([ 
    13.938 , 13.816,  13.661,  13.683,
    13.621,  13.547,  13.477,  13.492,
    13.237,  13.232,  13.07,   13.048,  
    12.945,  12.861,  12.827,  12.577,  
    12.518
])

y500 = np.array([  
    7.65103978e-06,   4.79865790e-06,   2.08218909e-05,   4.98385924e-06,
    5.63462673e-06,   2.90785458e-06,   2.21166794e-05,   1.34501705e-06,
    6.26021870e-07,   6.62368879e-07,   6.46735547e-07,   3.68589447e-07,
    3.86209019e-07,   5.61293275e-07,   2.41428755e-07,   9.62491134e-08,
    2.36892162e-07
])

代码

选项 1：绘制转换后的数据

# Regression Function
def regress(x, y):
    """Return a tuple of predicted y values and parameters for linear regression."""
    p = sp.stats.linregress(x, y)
    b1, b0, r, p_val, stderr = p
    y_pred = sp.polyval([b1, b0], x)
    return y_pred, p

# Plotting
x, y = mass500, np.log(y500)                      # transformed data
y_pred, _ = regress(x, y)

plt.plot(x, y, "mo", label="Data")
plt.plot(x, y_pred, "k--", label="Pred.")
plt.xlabel("Mass500")
plt.ylabel("log y500")                            # label axis
plt.legend()

输出

一种简单的方法是绘制转换后的数据并标记适当的对数轴。

选项 2：重新缩放轴并转换输入/输出函数值

代码

x, y = mass500, y500                              # data, non-transformed
y_pred, _ = regress(x, np.log(y))                 # transformed input             

plt.plot(x, y, "o", label="Data")
plt.plot(x, np.exp(y_pred), "k--", label="Pred.") # transformed output
plt.xlabel("Mass500")
plt.ylabel("y500")
plt.semilogy()
plt.legend()

输出

第二个选项是将轴更改为半对数刻度（通过plt.semilogy()）。在这里，未转换的数据自然呈现线性。另请注意，标签按原样表示数据。

要进行准确的回归，剩下的就是转换传递给回归函数的数据（通过np.log(x) 或np.log10(x)）以便返回正确的回归参数。当使用互补运算（即np.exp(x) 或10**x）绘制预测值时，这种转换会立即反转。

【讨论】：

【解决方案2】：

如果您想要一条在 log-y 尺度上看起来不错的线，则将一条线拟合到 y 值的对数。

def line(x, a, b):
    return a*x+b
popt, pcov = curve_fit(line, Mass500, np.log10(y500))
plt.semilogy(Mass500, 10**line(Mass500, popt[0], popt[1]), 'r-')

就是这样；我只省略了看起来不相关的三次插值部分。

【讨论】：