python - Pythonで指数曲線と対数曲線のフィッティングを行う方法は？多項式フィッティングしか見つかりませんでした

Question

私は一連のデータを持っており、どの行がそれを最もよく表しているかを比較したいと思います（異なる次数の多項式、指数関数または対数）。

私はPythonとNumpyを使用しており、多項式フィッティングには関数がありますpolyfit()。しかし、指数フィッティングと対数フィッティングのためのそのような関数は見つかりませんでした。

いずれかがあります？またはそれ以外の方法でそれを解決する方法は？

Answer 1

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

---

y = Ae ^Bxをフィッティングする場合、両側の対数をとると、log y = log A + Bxが得られます。したがって、 xに対して（log y）を当てはめます。

線形であるかのようにフィッティング（log y ）すると、 yの小さな値が強調され、大きなyに対して大きな偏差が発生することに注意してください。これは、polyfit（線形回帰）が∑ <sub> i（ΔY）² = ∑ <sub> i（Y i −Ŷi）2を最小_化する_こと^{によって機能するため}です。Y _i = log y _iの場合、剰余ΔY _i =Δ（logy _i）≈Δyy _/ | yi _|。だからpolyfit大きなy、「divide-by- | y |」に対して、非常に悪い決定を下します。ファクターはそれを補正し、polyfit小さな値を優先します。

これは、各エントリにyに比例する「重み」を与えることで軽減できます。キーワード引数polyfitを介して加重最小二乗法をサポートします。w

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

Excel、LibreOffice、およびほとんどの関数電卓は、通常、指数回帰/傾向線に重み付けされていない（バイアスされた）式を使用することに注意してください。結果をこれらのプラットフォームと互換性があるようにする場合は、より良い結果が得られる場合でも、重みを含めないでください。

---

これで、scipyを使用できる場合は、scipy.optimize.curve_fit変換なしで任意のモデルに適合させるために使用できます。

y = A + B log xの場合、結果は変換方法と同じです。

>>> x = numpy.array([1, 7, 20, 50, 79])
>>> y = numpy.array([10, 19, 30, 35, 51])
>>> scipy.optimize.curve_fit(lambda t,a,b: a+b*numpy.log(t),  x,  y)
(array([ 6.61867467,  8.46295606]), 
 array([[ 28.15948002,  -7.89609542],
        [ -7.89609542,   2.9857172 ]]))
# y ≈ 6.62 + 8.46 log(x)

ただし、 y = Ae ^Bxの場合、Δ（log y）を直接計算するため、より適切に適合させることができます。curve_fitただし、目的の極小値に到達できるように、初期化の推測を提供する必要があります。

>>> x = numpy.array([10, 19, 30, 35, 51])
>>> y = numpy.array([1, 7, 20, 50, 79])
>>> scipy.optimize.curve_fit(lambda t,a,b: a*numpy.exp(b*t),  x,  y)
(array([  5.60728326e-21,   9.99993501e-01]),
 array([[  4.14809412e-27,  -1.45078961e-08],
        [ -1.45078961e-08,   5.07411462e+10]]))
# oops, definitely wrong.
>>> scipy.optimize.curve_fit(lambda t,a,b: a*numpy.exp(b*t),  x,  y,  p0=(4, 0.1))
(array([ 4.88003249,  0.05531256]),
 array([[  1.01261314e+01,  -4.31940132e-02],
        [ -4.31940132e-02,   1.91188656e-04]]))
# y ≈ 4.88 exp(0.0553 x). much better.

Answer 2

curve_fitから使用する任意の関数にデータのセットを適合させることもできますscipy.optimize。たとえば、（ドキュメントからの）指数関数を適合させたい場合：

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)

そして、プロットしたい場合は、次のことができます。

plt.figure()
plt.plot(x, yn, 'ko', label="Original Noised Data")
plt.plot(x, func(x, *popt), 'r-', label="Fitted Curve")
plt.legend()
plt.show()

（注：プロットするとき*の前に、用語が、、に展開されます。これは予期されていることです。）poptabcfunc

Answer 3

私はこれに問題を抱えていたので、私のような初心者が理解できるように、非常に明確にさせてください。

データファイルなどがあるとしましょう

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

結果は次のとおりです。a=0.849195983017、b = -1.18101681765、c = 2.24061176543、d = 0.816643894816

Answer 4

これは、 scikitlearnのツールを使用する単純なデータの線形化オプションです。

与えられた

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

transformer = FunctionTransformer(np.log, validate=True)

コード

指数データを近似する

# Data
x_samp, y_samp = generate_data(func_exp, 2.5, 1.2, 0.7, jitter=3)
y_trans = transformer.fit_transform(y_samp)             # 1

# Regression
regressor = LinearRegression()
results = regressor.fit(x_samp, y_trans)                # 2
model = results.predict
y_fit = model(x_samp)

# Visualization
plt.scatter(x_samp, y_samp)
plt.plot(x_samp, np.exp(y_fit), "k--", label="Fit")     # 3
plt.title("Exponential Fit")

ログデータを適合させる

# Data
x_samp, y_samp = generate_data(func_log, 2.5, 1.2, 0.7, jitter=0.15)
x_trans = transformer.fit_transform(x_samp)             # 1

# Regression
regressor = LinearRegression()
results = regressor.fit(x_trans, y_samp)                # 2
model = results.predict
y_fit = model(x_trans)

# Visualization
plt.scatter(x_samp, y_samp)
plt.plot(x_samp, y_fit, "k--", label="Fit")             # 3
plt.title("Logarithmic Fit")

---

詳細

一般的な手順

1. データ値（x、yまたは両方）にログ操作を適用します
2. データを線形化されたモデルに回帰します
3. （を使用して）ログ操作を「逆」にしてプロットしnp.exp()、元のデータに適合させます

データが指数関数的な傾向に従うと仮定すると、一般的な方程式⁺は次のようになります。

後者の方程式（たとえば、y=切片+傾き*x）は、対数を取ることで線形化できます。

^{線形化された方程式++}と回帰パラメーターが与えられると、次のように計算できます。

• Aインターセプト経由（ln(A)）
• Bスロープ経由（B）

線形化手法の概要

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)

_{⁺注：指数関数の線形化は、ノイズが小さく、C=0の場合に最適に機能します。注意して使用してください。}

_{⁺⁺注：xデータを変更すると指数データを線形化できますが、yデータを変更するとログデータを線形化できます。}

Answer 5

さて、私はあなたがいつでも使うことができると思います：

np.log   -->  natural log
np.log10 -->  base 10
np.log2  -->  base 2

---

IanVSの答えを少し変更します：

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
  return a * np.log(b * x) + c

x = np.linspace(1,5,50)   # changed boundary conditions to avoid division by 0
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)

plt.figure()
plt.plot(x, yn, 'ko', label="Original Noised Data")
plt.plot(x, func(x, *popt), 'r-', label="Fitted Curve")
plt.legend()
plt.show()

これにより、次のグラフが作成されます。

Answer 6

lmfit両方の問題を解決しながらの機能を示します。

与えられた

import lmfit

import numpy as np

import matplotlib.pyplot as plt


%matplotlib inline
np.random.seed(123)

# General Functions
def 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

コード

アプローチ1-lmfitモデル

指数データを近似する

regressor = lmfit.models.ExponentialModel()                # 1    
initial_guess = dict(amplitude=1, decay=-1)                # 2
results = regressor.fit(y_samp, x=x_samp, **initial_guess)
y_fit = results.best_fit    

plt.plot(x_samp, y_samp, "o", label="Data")
plt.plot(x_samp, y_fit, "k--", label="Fit")
plt.legend()

アプローチ2-カスタムモデル

ログデータを適合させる

regressor = lmfit.Model(func_log)                          # 1
initial_guess = dict(a=1, b=.1, c=.1)                      # 2
results = regressor.fit(y_samp2, x=x_samp, **initial_guess)
y_fit = results.best_fit

plt.plot(x_samp, y_samp2, "o", label="Data")
plt.plot(x_samp, y_fit, "k--", label="Fit")
plt.legend()

---

詳細

1. 回帰クラスを選択する
2. 関数の定義域を尊重する名前付きの初期推測を提供します

リグレッサオブジェクトから推測されたパラメータを決定できます。例：

regressor.param_names
# ['decay', 'amplitude']

予測を行うには、このメソッドを使用しModelResult.eval()ます。

model = results.eval
y_pred = model(x=np.array([1.5]))

注：ExponentialModel()以下は、2つのパラメーターを受け入れる減衰関数であり、そのうちの1つは負です。

より多くのパラメータExponentialGaussianModel()を受け入れる、も参照してください。

を介してライブラリをインストール> pip install lmfitします。

Answer 7

Wolframには、指数関数を近似するための閉じた形の解があります。また、対数法則とべき法則を適合させるための同様のソリューションもあります。

これはscipyのcurve_fitよりもうまく機能することがわかりました。特に「ゼロに近い」データがない場合。次に例を示します。

import numpy as np
import matplotlib.pyplot as plt

# Fit the function y = A * exp(B * x) to the data
# returns (A, B)
# From: https://mathworld.wolfram.com/LeastSquaresFittingExponential.html
def fit_exp(xs, ys):
    S_x2_y = 0.0
    S_y_lny = 0.0
    S_x_y = 0.0
    S_x_y_lny = 0.0
    S_y = 0.0
    for (x,y) in zip(xs, ys):
        S_x2_y += x * x * y
        S_y_lny += y * np.log(y)
        S_x_y += x * y
        S_x_y_lny += x * y * np.log(y)
        S_y += y
    #end
    a = (S_x2_y * S_y_lny - S_x_y * S_x_y_lny) / (S_y * S_x2_y - S_x_y * S_x_y)
    b = (S_y * S_x_y_lny - S_x_y * S_y_lny) / (S_y * S_x2_y - S_x_y * S_x_y)
    return (np.exp(a), b)


xs = [33, 34, 35, 36, 37, 38, 39, 40, 41, 42]
ys = [3187, 3545, 4045, 4447, 4872, 5660, 5983, 6254, 6681, 7206]

(A, B) = fit_exp(xs, ys)

plt.figure()
plt.plot(xs, ys, 'o-', label='Raw Data')
plt.plot(xs, [A * np.exp(B *x) for x in xs], 'o-', label='Fit')

plt.title('Exponential Fit Test')
plt.xlabel('X')
plt.ylabel('Y')
plt.legend(loc='best')
plt.tight_layout()
plt.show()