I have a set of data and I want to compare which line describes it best (polynomials of different orders, exponential or logarithmic).
I use Python and Numpy and for polynomial fitting there is a function polyfit(). But I found no such functions for exponential and logarithmic fitting.
Are there any? Or how to solve it otherwise?
回答 0
对于拟合y = A + B log x,只需将y拟合为(log x)。
>>> x = numpy.array([1,7,20,50,79])>>> y = numpy.array([10,19,30,35,51])>>> numpy.polyfit(numpy.log(x), y,1)
array([8.46295607,6.61867463])# y ≈ 8.46 log(x) + 6.62
For fitting y = A + B log x, just fit y against (log x).
>>> x = numpy.array([1, 7, 20, 50, 79])
>>> y = numpy.array([10, 19, 30, 35, 51])
>>> numpy.polyfit(numpy.log(x), y, 1)
array([ 8.46295607, 6.61867463])
# y ≈ 8.46 log(x) + 6.62
For fitting y = AeBx, take the logarithm of both side gives log y = log A + Bx. So fit (log y) against x.
Note that fitting (log y) as if it is linear will emphasize small values of y, causing large deviation for large y. This is because polyfit (linear regression) works by minimizing ∑i (ΔY)2 = ∑i (Yi − Ŷi)2. When Yi = log yi, the residues ΔYi = Δ(log yi) ≈ Δyi / |yi|. So even if polyfit makes a very bad decision for large y, the “divide-by-|y|” factor will compensate for it, causing polyfit favors small values.
This could be alleviated by giving each entry a “weight” proportional to y. polyfit supports weighted-least-squares via the w keyword argument.
>>> x = numpy.array([10, 19, 30, 35, 51])
>>> y = numpy.array([1, 7, 20, 50, 79])
>>> numpy.polyfit(x, numpy.log(y), 1)
array([ 0.10502711, -0.40116352])
# y ≈ exp(-0.401) * exp(0.105 * x) = 0.670 * exp(0.105 * x)
# (^ biased towards small values)
>>> numpy.polyfit(x, numpy.log(y), 1, w=numpy.sqrt(y))
array([ 0.06009446, 1.41648096])
# y ≈ exp(1.42) * exp(0.0601 * x) = 4.12 * exp(0.0601 * x)
# (^ not so biased)
Note that Excel, LibreOffice and most scientific calculators typically use the unweighted (biased) formula for the exponential regression / trend lines. If you want your results to be compatible with these platforms, do not include the weights even if it provides better results.
Now, if you can use scipy, you could use scipy.optimize.curve_fit to fit any model without transformations.
For y = A + B log x the result is the same as the transformation method:
For y = AeBx, however, we can get a better fit since it computes Δ(log y) directly. But we need to provide an initialize guess so curve_fit can reach the desired local minimum.
import numpy as npimport matplotlib.pyplot as pltfrom scipy.optimize import curve_fitdef func(x, a, b, c):return a * np.exp(-b * x)+ c
x = np.linspace(0,4,50)
y = func(x,2.5,1.3,0.5)
yn = y +0.2*np.random.normal(size=len(x))
popt, pcov = curve_fit(func, x, yn)
You can also fit a set of a data to whatever function you like using curve_fit from scipy.optimize. For example if you want to fit an exponential function (from the documentation):
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
def func(x, a, b, c):
return a * np.exp(-b * x) + c
x = np.linspace(0,4,50)
y = func(x, 2.5, 1.3, 0.5)
yn = y + 0.2*np.random.normal(size=len(x))
popt, pcov = curve_fit(func, x, yn)
(Note: the * in front of popt when you plot will expand out the terms into the a, b, and c that func is expecting.)
回答 2
我对此有些麻烦,所以请让我非常明确,让像我这样的菜鸟可以理解。
假设我们有一个数据文件或类似的文件
# -*- coding: utf-8 -*-import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import numpy as np
import sympy as sym
"""
Generate some data, let's imagine that you already have this.
"""
x = np.linspace(0,3,50)
y = np.exp(x)"""
Plot your data
"""
plt.plot(x, y,'ro',label="Original Data")"""
brutal force to avoid errors
"""
x = np.array(x, dtype=float)#transform your data in a numpy array of floats
y = np.array(y, dtype=float)#so the curve_fit can work"""
create a function to fit with your data. a, b, c and d are the coefficients
that curve_fit will calculate for you.
In this part you need to guess and/or use mathematical knowledge to find
a function that resembles your data
"""def func(x, a, b, c, d):return a*x**3+ b*x**2+c*x + d
"""
make the curve_fit
"""
popt, pcov = curve_fit(func, x, y)"""
The result is:
popt[0] = a , popt[1] = b, popt[2] = c and popt[3] = d of the function,
so f(x) = popt[0]*x**3 + popt[1]*x**2 + popt[2]*x + popt[3].
"""print"a = %s , b = %s, c = %s, d = %s"%(popt[0], popt[1], popt[2], popt[3])"""
Use sympy to generate the latex sintax of the function
"""
xs = sym.Symbol('\lambda')
tex = sym.latex(func(xs,*popt)).replace('$','')
plt.title(r'$f(\lambda)= %s$'%(tex),fontsize=16)"""
Print the coefficients and plot the funcion.
"""
plt.plot(x, func(x,*popt), label="Fitted Curve")#same as line above \/#plt.plot(x, popt[0]*x**3 + popt[1]*x**2 + popt[2]*x + popt[3], label="Fitted Curve")
plt.legend(loc='upper left')
plt.show()
I was having some trouble with this so let me be very explicit so noobs like me can understand.
Lets say that we have a data file or something like that
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import numpy as np
import sympy as sym
"""
Generate some data, let's imagine that you already have this.
"""
x = np.linspace(0, 3, 50)
y = np.exp(x)
"""
Plot your data
"""
plt.plot(x, y, 'ro',label="Original Data")
"""
brutal force to avoid errors
"""
x = np.array(x, dtype=float) #transform your data in a numpy array of floats
y = np.array(y, dtype=float) #so the curve_fit can work
"""
create a function to fit with your data. a, b, c and d are the coefficients
that curve_fit will calculate for you.
In this part you need to guess and/or use mathematical knowledge to find
a function that resembles your data
"""
def func(x, a, b, c, d):
return a*x**3 + b*x**2 +c*x + d
"""
make the curve_fit
"""
popt, pcov = curve_fit(func, x, y)
"""
The result is:
popt[0] = a , popt[1] = b, popt[2] = c and popt[3] = d of the function,
so f(x) = popt[0]*x**3 + popt[1]*x**2 + popt[2]*x + popt[3].
"""
print "a = %s , b = %s, c = %s, d = %s" % (popt[0], popt[1], popt[2], popt[3])
"""
Use sympy to generate the latex sintax of the function
"""
xs = sym.Symbol('\lambda')
tex = sym.latex(func(xs,*popt)).replace('$', '')
plt.title(r'$f(\lambda)= %s$' %(tex),fontsize=16)
"""
Print the coefficients and plot the funcion.
"""
plt.plot(x, func(x, *popt), label="Fitted Curve") #same as line above \/
#plt.plot(x, popt[0]*x**3 + popt[1]*x**2 + popt[2]*x + popt[3], label="Fitted Curve")
plt.legend(loc='upper left')
plt.show()
the result is:
a = 0.849195983017 , b = -1.18101681765, c = 2.24061176543, d = 0.816643894816
回答 3
好吧,我想您可以随时使用:
np.log --> natural log
np.log10 --> base 10
np.log2 --> base 2
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model importLinearRegressionfrom sklearn.preprocessing importFunctionTransformer
np.random.seed(123)
# General Functionsdef func_exp(x, a, b, c):"""Return values from a general exponential function."""return a * np.exp(b * x)+ c
def func_log(x, a, b, c):"""Return values from a general log function."""return a * np.log(b * x)+ c
# Helperdef generate_data(func,*args, jitter=0):"""Return a tuple of arrays with random data along a general function."""
xs = np.linspace(1,5,50)
ys = func(xs,*args)
noise = jitter * np.random.normal(size=len(xs))+ jitter
xs = xs.reshape(-1,1)# xs[:, np.newaxis]
ys =(ys + noise).reshape(-1,1)return xs, ys
Relationship|Example|GeneralEqn.|AlteredVar.|LinearizedEqn.-------------|------------|----------------------|----------------|------------------------------------------Linear| x | y = B * x + C |-| y = C + B * x
Logarithmic| log(x)| y = A * log(B*x)+ C | log(x)| y = C + A *(log(B)+ log(x))Exponential|2**x, e**x | y = A * exp(B*x)+ C | log(y)| log(y-C)= log(A)+ B * x
Power| x**2| y = B * x**N + C | log(x), log(y)| log(y-C)= log(B)+ N * log(x)
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import FunctionTransformer
np.random.seed(123)
# General Functions
def func_exp(x, a, b, c):
"""Return values from a general exponential function."""
return a * np.exp(b * x) + c
def func_log(x, a, b, c):
"""Return values from a general log function."""
return a * np.log(b * x) + c
# Helper
def generate_data(func, *args, jitter=0):
"""Return a tuple of arrays with random data along a general function."""
xs = np.linspace(1, 5, 50)
ys = func(xs, *args)
noise = jitter * np.random.normal(size=len(xs)) + jitter
xs = xs.reshape(-1, 1) # xs[:, np.newaxis]
ys = (ys + noise).reshape(-1, 1)
return xs, ys
Apply a log operation to data values (x, y or both)
Regress the data to a linearized model
Plot by “reversing” any log operations (with np.exp()) and fit to original data
Assuming our data follows an exponential trend, a general equation+ may be:
We can linearize the latter equation (e.g. y = intercept + slope * x) by taking the log:
Given a linearized equation++ and the regression parameters, we could calculate:
A via intercept (ln(A))
B via slope (B)
Summary of Linearization Techniques
Relationship | Example | General Eqn. | Altered Var. | Linearized Eqn.
-------------|------------|----------------------|----------------|------------------------------------------
Linear | x | y = B * x + C | - | y = C + B * x
Logarithmic | log(x) | y = A * log(B*x) + C | log(x) | y = C + A * (log(B) + log(x))
Exponential | 2**x, e**x | y = A * exp(B*x) + C | log(y) | log(y-C) = log(A) + B * x
Power | x**2 | y = B * x**N + C | log(x), log(y) | log(y-C) = log(B) + N * log(x)
+Note: linearizing exponential functions works best when the noise is small and C=0. Use with caution.
++Note: while altering x data helps linearize exponential data, altering y data helps linearize log data.
import lmfit
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed(123)
# General Functionsdef func_log(x, a, b, c):"""Return values from a general log function."""return a * np.log(b * x)+ c
# Data
x_samp = np.linspace(1,5,50)
_noise = np.random.normal(size=len(x_samp), scale=0.06)
y_samp =2.5* np.exp(1.2* x_samp)+0.7+ _noise
y_samp2 =2.5* np.log(1.2* x_samp)+0.7+ _noise