拟合曲线以获取由两个不同的方案组成的数据 [英] Fit a curve for data made up of two distinct regimes

问题描述

我正在寻找一种通过一些实验数据绘制曲线的方法.数据显示出一个小的线性方案，其梯度较浅，随后在阈值之后为陡峭的线性方案.

我的数据在这里: http://pastebin.com/H4NSbxqr

我可以相对容易地用两条线拟合数据，但我想理想地用连续线拟合-看起来应该像两条平滑的曲线将它们连接在阈值附近(在数据中为〜5000，如图所示)以上).

我尝试使用scipy.optimize curve_fit并尝试包含直线和指数之和的函数:

y = a*x + b + c*np.exp((x-d)/e)

尽管进行了许多尝试，但仍未找到解决方案.

如果有人对拟合分布/方法的选择或curve_fit实现有任何建议，将不胜感激.

解决方案

如果您没有特殊的理由认为线性+指数是数据的真正根本原因，那么我认为适合两行最有意义的.您可以通过使拟合函数最大为两行来实现此目的，例如:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def two_lines(x, a, b, c, d):
    one = a*x + b
    two = c*x + d
    return np.maximum(one, two)

然后

x, y = np.genfromtxt('tmp.txt', unpack=True, delimiter=',')

pw0 = (.02, 30, .2, -2000) # a guess for slope, intercept, slope, intercept
pw, cov = curve_fit(two_lines, x, y, pw0)
crossover = (pw[3] - pw[1]) / (pw[0] - pw[2])

plt.plot(x, y, 'o', x, two_lines(x, *pw), '-')

如果您真的想要一个连续且有区别的解决方案，对我来说，双曲线有一个急剧的弯曲，但它必须旋转.实施起来有点困难(也许有一种更简单的方法)，但是可以尝试一下:

def hyperbola(x, a, b, c, d, e):
    """ hyperbola(x) with parameters
        a/b = asymptotic slope
         c  = curvature at vertex
         d  = offset to vertex
         e  = vertical offset
    """
    return a*np.sqrt((b*c)**2 + (x-d)**2)/b + e

def rot_hyperbola(x, a, b, c, d, e, th):
    pars = a, b, c, 0, 0 # do the shifting after rotation
    xd = x - d
    hsin = hyperbola(xd, *pars)*np.sin(th)
    xcos = xd*np.cos(th)
    return e + hyperbola(xcos - hsin, *pars)*np.cos(th) + xcos - hsin

运行为

h0 = 1.1, 1, 0, 5000, 100, .5
h, hcov = curve_fit(rot_hyperbola, x, y, h0)
plt.plot(x, y, 'o', x, two_lines(x, *pw), '-', x, rot_hyperbola(x, *h), '-')
plt.legend(['data', 'piecewise linear', 'rotated hyperbola'], loc='upper left')
plt.show()

我也能够使线+指数收敛，但是看起来很糟糕.这是因为它不是很好的数据描述符，它是线性的，而指数则远远不是线性的！

def line_exp(x, a, b, c, d, e):
    return a*x + b + c*np.exp((x-d)/e)

e0 = .1, 20., .01, 1000., 2000.
e, ecov = curve_fit(line_exp, x, y, e0)

如果要保持简单，总会有多项式或样条曲线(分段多项式)

from scipy.interpolate import UnivariateSpline
s = UnivariateSpline(x, y, s=x.size)  #larger s-value has fewer "knots"
plt.plot(x, s(x))

I'm looking for a way to plot a curve through some experimental data. The data shows a small linear regime with a shallow gradient, followed by a steep linear regime after a threshold value.

My data is here: http://pastebin.com/H4NSbxqr

I could fit the data with two lines relatively easily, but I'd like to fit with a continuous line ideally - which should look like two lines with a smooth curve joining them around the threshold (~5000 in the data, shown above).

I attempted this using scipy.optimize curve_fit and trying a function which included the sum of a straight line and an exponential:

y = a*x + b + c*np.exp((x-d)/e)

although despite numerous attempts, it didn't find a solution.

If anyone has any suggestions please, either on the choice of fitting distribution / method or the curve_fit implementation, they would be greatly appreciated.

解决方案

If you don't have a particular reason to believe that linear + exponential is the true underlying cause of your data, then I think a fit to two lines makes the most sense. You can do this by making your fitting function the maximum of two lines, for example:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def two_lines(x, a, b, c, d):
    one = a*x + b
    two = c*x + d
    return np.maximum(one, two)

Then,

x, y = np.genfromtxt('tmp.txt', unpack=True, delimiter=',')

pw0 = (.02, 30, .2, -2000) # a guess for slope, intercept, slope, intercept
pw, cov = curve_fit(two_lines, x, y, pw0)
crossover = (pw[3] - pw[1]) / (pw[0] - pw[2])

plt.plot(x, y, 'o', x, two_lines(x, *pw), '-')

If you really want a continuous and differentiable solution, it occurred to me that a hyperbola has a sharp bend to it, but it has to be rotated. It was a bit difficult to implement (maybe there's an easier way), but here's a go:

def hyperbola(x, a, b, c, d, e):
    """ hyperbola(x) with parameters
        a/b = asymptotic slope
         c  = curvature at vertex
         d  = offset to vertex
         e  = vertical offset
    """
    return a*np.sqrt((b*c)**2 + (x-d)**2)/b + e

def rot_hyperbola(x, a, b, c, d, e, th):
    pars = a, b, c, 0, 0 # do the shifting after rotation
    xd = x - d
    hsin = hyperbola(xd, *pars)*np.sin(th)
    xcos = xd*np.cos(th)
    return e + hyperbola(xcos - hsin, *pars)*np.cos(th) + xcos - hsin

Run it as

h0 = 1.1, 1, 0, 5000, 100, .5
h, hcov = curve_fit(rot_hyperbola, x, y, h0)
plt.plot(x, y, 'o', x, two_lines(x, *pw), '-', x, rot_hyperbola(x, *h), '-')
plt.legend(['data', 'piecewise linear', 'rotated hyperbola'], loc='upper left')
plt.show()

I was also able to get the line + exponential to converge, but it looks terrible. This is because it's not a good descriptor of your data, which is linear and an exponential is very far from linear!

def line_exp(x, a, b, c, d, e):
    return a*x + b + c*np.exp((x-d)/e)

e0 = .1, 20., .01, 1000., 2000.
e, ecov = curve_fit(line_exp, x, y, e0)

If you want to keep it simple, there's always a polynomial or spline (piecewise polynomials)

from scipy.interpolate import UnivariateSpline
s = UnivariateSpline(x, y, s=x.size)  #larger s-value has fewer "knots"
plt.plot(x, s(x))

这篇关于拟合曲线以获取由两个不同的方案组成的数据的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！