Cara menggunakan SCIPY.OPTIMIZE pada Python

Saya memiliki satu set data dan saya ingin membandingkan baris mana yang paling menggambarkannya (polinomial dengan urutan berbeda, eksponensial, atau logaritmik).

Saya menggunakan Python dan Numpy dan untuk pemasangan polinom ada fungsi polyfit(). Tetapi saya tidak menemukan fungsi seperti itu untuk pemasangan eksponensial dan logaritmik.

Apakah ada? Atau bagaimana cara menyelesaikannya?

Untuk pemasangan y = A + B log x , cukup pas y terhadap (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

Untuk pemasangan y = Ae Bx , ambil logaritma dari kedua sisi beri log y = log A + Bx . Jadi pas (log y ) terhadap x .

Perhatikan bahwa pemasangan (log y ) seolah-olah linier akan menekankan nilai-nilai kecil y , menyebabkan penyimpangan besar untuk y besar . Ini karena polyfit(regresi linier) bekerja dengan meminimalkan ∑ i (Δ Y ) 2 = ∑ i ( Y i - Ŷ i ) 2 . Ketika Y i = log y i , residu Δ Y i = Δ (log y i ) ≈ Δ y i / | y i |. Bahkan jikapolyfitmembuat keputusan yang sangat buruk untuk y besar , "divide-by- | y |" Faktor akan mengimbanginya, menyebabkan polyfitnilai-nilai kecil nikmat.

Ini dapat dikurangi dengan memberikan setiap entri "bobot" sebanding dengan y . polyfitmendukung weighted-least-square melalui wargumen kata kunci.

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

Perhatikan bahwa Excel, LibreOffice, dan sebagian besar kalkulator ilmiah biasanya menggunakan rumus tanpa bobot (bias) untuk garis regresi / tren eksponensial. Jika Anda ingin hasil Anda kompatibel dengan platform ini, jangan masukkan bobot meskipun itu memberikan hasil yang lebih baik.

Sekarang, jika Anda bisa menggunakan Scipy, Anda bisa menggunakan scipy.optimize.curve_fitmodel apa saja tanpa transformasi.

Untuk y = log A + B x hasilnya sama dengan metode transformasi:

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

Namun untuk y = Ae Bx , kita bisa mendapatkan kecocokan yang lebih baik karena ia menghitung Δ (log y ) secara langsung. Tetapi kita perlu memberikan perkiraan inisialisasi sehingga curve_fitdapat mencapai minimum lokal yang diinginkan.

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

Anda juga dapat muat satu set data untuk apa pun yang berfungsi Anda suka menggunakan curve_fitdari scipy.optimize. Misalnya jika Anda ingin menyesuaikan fungsi eksponensial (dari dokumentasi ):

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)

Dan jika Anda ingin merencanakan, Anda bisa melakukan:

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

(Catatan: *di depan poptketika Anda merencanakan akan memperluas syarat-syarat ke dalam a, bdan cyang func. Mengharapkan)

Saya mengalami beberapa masalah dengan ini jadi biarkan saya menjadi sangat eksplisit sehingga noobs seperti saya bisa mengerti.

Katakanlah kita memiliki file data atau sesuatu seperti itu

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

hasilnya adalah: a = 0.849195983017, b = -1.18101681765, c = 2.24061176543, d = 0.816643894816

Yah saya kira Anda selalu dapat menggunakan:

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

Sedikit memodifikasi jawaban 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()

Ini menghasilkan grafik berikut:

Berikut adalah opsi linierisasi pada data sederhana yang menggunakan alat dari scikit belajar .

Diberikan

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)

Kode

Sesuai dengan data eksponensial

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

Sesuaikan data log

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

Detail

Langkah Umum

1. Terapkan operasi log ke nilai data ( x, yatau keduanya)
2. Regress data ke model linier
3. Plot dengan "membalikkan" operasi log apa pun (dengan np.exp()) dan pas dengan data asli

Dengan asumsi data kami mengikuti tren eksponensial, persamaan umum + mungkin:

Kita dapat meratakan persamaan kedua (misalnya y = intersep + kemiringan * x) dengan mengambil log :

Diberikan persamaan linear + dan parameter regresi, kita dapat menghitung:

• Avia intercept ( ln(A))
• Bvia slope ( B)

Ringkasan Teknik Linearisasi

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)

+ Catatan: fungsi eksponensial linier bekerja paling baik saat noise kecil dan C = 0. Gunakan dengan hati-hati.

++ Catatan: sementara mengubah data x membantu linierisasi data eksponensial , mengubah data y membantu linierisasi data log .

Kami menunjukkan fitur lmfitsambil menyelesaikan kedua masalah.

Diberikan

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

Kode

Pendekatan 1 - lmfitModel

Sesuai dengan data eksponensial

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

Pendekatan 2 - Model Kustom

Sesuaikan data log

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

Detail

1. Pilih kelas regresi
2. Nama pasokan, tebakan awal yang menghormati domain fungsi

Anda dapat menentukan parameter yang disimpulkan dari objek regressor. Contoh:

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

Catatan: ExponentialModel()berikut ini adalah fungsi peluruhan , yang menerima dua parameter, salah satunya negatif.

Lihat juga ExponentialGaussianModel(), yang menerima lebih banyak parameter .

Instal perpustakaan melalui > pip install lmfit.

Wolfram memiliki solusi bentuk tertutup untuk memasang eksponensial . Mereka juga memiliki solusi serupa untuk menyesuaikan hukum logaritmik dan kekuasaan .

Saya menemukan ini bekerja lebih baik daripada curve_fit scipy's. Berikut ini sebuah contoh:

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