Hồi quy phi tuyến tính Python

Question

Kumar, A. , & Babcock, J. (2017). con trăn. Phân tích dự đoán nâng cao. Có được những hiểu biết thực tế bằng cách khai thác dữ liệu trong doanh nghiệp của bạn để xây dựng các ứng dụng mô hình dự đoán nâng cao. Công ty TNHH xuất bản Packt

Trong thống kê, chúng tôi nói rằng hồi quy là tuyến tính khi nó tuyến tính trong các tham số. Khớp các mô hình tuyến tính là một nhiệm vụ dễ dàng, chúng ta có thể sử dụng phương pháp bình phương nhỏ nhất và lấy các tham số tối ưu cho mô hình của mình. Trong Python, bạn có thể đạt được điều này bằng cách sử dụng một loạt các thư viện như

from sympy import *
a, b, c, x = symbols("a b c x")
expr = a*x**2 + b*x + c
diff(expr, a), diff(expr, b), diff(expr, c)

0,

from sympy import *
a, b, c, x = symbols("a b c x")
expr = a*x**2 + b*x + c
diff(expr, a), diff(expr, b), diff(expr, c)

1,

from sympy import *
a, b, c, x = symbols("a b c x")
expr = a*x**2 + b*x + c
diff(expr, a), diff(expr, b), diff(expr, c)

2,

from sympy import *
a, b, c, x = symbols("a b c x")
expr = a*x**2 + b*x + c
diff(expr, a), diff(expr, b), diff(expr, c)

3, v.v. Tuy nhiên, không phải tất cả các vấn đề có thể được giải quyết với các mô hình tuyến tính thuần túy. Có rất nhiều mô hình phi tuyến hữu ích đảm bảo các tính chất toán học hữu ích và cũng rất dễ hiểu. Vì tôi có xu hướng thường xuyên điều chỉnh các mô hình này, nên bài viết t sẽ đóng vai trò là tài liệu tham khảo trong tương lai cho tôi bất cứ khi nào tôi cần sử dụng chúng trong tương lai, với mã dễ triển khai (bằng Python)

Nội dung chính Show

Tuyến tính x Phi tuyến
Lắp đường cong trong Python
Đoán ban đầu và Jacobian
Mô hình lồi/lõm
Phân rã theo cấp số nhân
Phân rã theo cấp số nhân với tiệm cận thấp hơn
Mô hình tiệm cận (Số mũ âm)
Mô hình tiệm cận (bị ràng buộc. bắt đầu từ 0)
hồi quy sức mạnh
Đường cong sigmoidal
đường cong hậu cần
Chức năng Gompertz
Kết luận + Mã
không phải là gì
Python hồi quy tuyến tính và phi tuyến tính là gì?
không phải là gì
không phải là gì

Bây giờ, trước khi chúng ta đi vào các mô hình, điều quan trọng cần lưu ý là bối cảnh về lý do và khi nào bạn sẽ sử dụng các loại mô hình này. Chúng không nhất thiết phải thay thế cho tổ hợp của bạn (XGBoost, LightGBM, Random Forest, Neural Network, v.v. ) mô hình học máy. Tốt nhất là nghĩ về các mô hình phi tuyến tính này trong bối cảnh khớp đường cong, tôi. e. bạn muốn mô hình hóa một hiện tượng bằng cách sử dụng một đường cong mà sau này bạn có thể giải thích. Vì vậy, đây là một số cân nhắc về lý do tại sao bạn sẽ sử dụng chúng

khả năng diễn giải. Nếu bạn lập mô hình một sự kiện bằng cách sử dụng các đường cong, thì kết quả của bạn sẽ dễ hiểu hơn nhiều, vì thông thường phương trình toán học có ý nghĩa đối với từng tham số. Quan trọng. điều này không có nghĩa là bạn mất đi sức mạnh của mô hình máy học phức tạp của mình, vì bạn vẫn có thể sử dụng mô hình máy học cùng với đường cong mà bạn đang điều chỉnh (e. g. phân đoạn dữ liệu của bạn bằng mô hình, sau đó điều chỉnh các đường cong có thể hiểu được trên từng phân đoạn). Điều này sẽ giúp ích không chỉ cho người xây dựng mô hình trong tương lai mà còn cho các bên liên quan trong kinh doanh. Sẽ dễ dàng hơn nhiều để giải thích mô hình của bạn đang làm gì nếu đầu ra có thể hiểu được
Thuộc tính toán học. Đôi khi, điều thực sự quan trọng là đảm bảo một số thuộc tính toán học trên đầu ra mô hình của bạn. Ví dụ: nếu bạn đang lập mô hình
```
from sympy import *
a, b, c, x = symbols("a b c x")
expr = a*x**2 + b*x + c
diff(expr, a), diff(expr, b), diff(expr, c)
```
4 của một số loài côn trùng, thì có thể dự đoán chiều cao của bạn sẽ tăng một cách đơn điệu khi tuổi tăng lên. Việc mô tả những hiện tượng này bằng các phương trình toán học sẽ cho bạn khả năng áp dụng ràng buộc này trong giai đoạn lập mô hình. Bạn vừa khớp với một hàm có thuộc tính đơn điệu này

Tuyến tính x Phi tuyến

Điều quan trọng là phải phân biệt rằng khi tôi nói phi tuyến tính, ý tôi là phi tuyến tính trong các tham số. Vì vậy, ví dụ, một phương trình bậc hai $y = ax^2 + bx + c$ là tuyến tính, bởi vì mô hình này là tuyến tính trong các tham số $a, b, c$. Mô hình có thể không tuyến tính trong $x$, nhưng nó vẫn có thể tuyến tính trong các tham số

Để hiểu rõ hơn về các mô hình tuyến tính và phi tuyến tính, hãy xem xét các ví dụ sau

$$ \begin{equation} y = \beta_0 + \beta_1 x \end{equation} $$

$$ \begin{equation} y = \beta_0 (1 + \beta_1)^x \end{equation} $$

$$ \begin{equation} y = \ \beta_0\cdot\sin\left(x^{\beta_1}\right)\ +\ \beta_2\cdot\cos\left(e^{x\beta_3}\right)

Phương trình $(1)$ là một dòng đơn giản và các tham số $\beta_0, \beta_1$ là tuyến tính trên $y$, vì vậy đây là một ví dụ về mô hình tuyến tính. Phương trình $(2)$ là cái mà chúng tôi gọi là công thức Tăng trưởng theo cấp số nhân, trong đó $\beta_0$ thường biểu thị điểm bắt đầu, $x$ là thước đo thời gian và $\beta_1$ (thường được gọi là $r$) biểu thị tốc độ tăng trưởng. Trong trường hợp này, chúng tôi muốn có được điểm xuất phát và tốc độ tăng trưởng;

$$ y = \beta_0 (1 + \beta_1)^x \\\
log(y) = log(\beta_0) + x\cdot log(1 + \beta_1) $$

Và sau đó gọi $y'=log(y)$, $\beta_0’ = log(\beta_0)$, $\beta_1'=log(1 + \beta_1)$, bạn có thể viết lại Tăng trưởng theo cấp số nhân thành

$$y’ = \beta_0’ + \beta_1’x$$

Và điều chỉnh OLS (Bình phương nhỏ nhất thông thường) bằng cách sử dụng công thức này, vì đây là mô hình tuyến tính (đây được gọi là mô hình log-tuyến tính). Cuối cùng, phương trình $(3)$ là một mô hình phi tuyến tính, vì hồi quy là phi tuyến tính trong các tham số $\beta_0, \beta_1, \beta_2, \beta_3, \beta_4$ (Tôi hy vọng vậy, khi tôi tạo nhanh phương trình này

Trong bài viết này, trước tiên tôi sẽ trình bày các công cụ chính mà chúng ta sẽ sử dụng để điều chỉnh các mô hình, sau đó giải thích một loạt các mô hình phi tuyến hữu ích + mã + đồ thị của mô hình cho bất cứ khi nào bạn cần điều chỉnh các phương trình kỳ diệu như vậy. Tôi cũng sẽ không đi sâu hơn vào mục đích và cách sử dụng từng model vì nó không phải trọng tâm của bài viết.

Lắp đường cong trong Python

Một phần lớn của bài viết này dựa trên một blogpost khác, nhằm vào người dùng R và sử dụng các thư viện R với các phương trình được tạo sẵn và các quy trình tích hợp sẵn để tạo điều kiện thuận lợi cho cuộc sống của nhà khoa học dữ liệu. Mặc dù đôi khi tôi sử dụng R, nhưng hầu hết thời gian tôi sử dụng Python, vì vậy thật tuyệt khi giữ cơ sở mã dự án của tôi với số lượng ngôn ngữ lập trình riêng biệt tối thiểu. Đáng buồn thay, không có thư viện dễ sử dụng/phổ biến nào trong Python tích hợp sẵn các mô hình này (nếu bạn biết bất kỳ thư viện nào, vui lòng cho tôi biết. ), nên chúng ta sẽ phải thực hiện các mô hình và sáng tạo trong quá trình điều chỉnh đường cong khi cần thiết, đặc biệt là khi chọn các tham số ban đầu.

Bên cạnh bộ ba thần thánh (Pandas, Matplotlib và Numpy), đây là những thư viện chính mà chúng tôi sẽ sử dụng, cho các mục đích phù hợp

scipy
Sympy

Chúng tôi chủ yếu sẽ sử dụng các quy trình tối ưu hóa từ scipy's

from sympy import *
a, b, c, x = symbols("a b c x")
expr = a*x**2 + b*x + c
diff(expr, a), diff(expr, b), diff(expr, c)

5. Các tham số chính chúng tôi sẽ sử dụng là

from sympy import *
a, b, c, x = symbols("a b c x")
expr = a*x**2 + b*x + c
diff(expr, a), diff(expr, b), diff(expr, c)

6. Chức năng gọi khi tối ưu hóa, tôi. e. mô hình của bạn. Nó phải có chữ ký tương tự như

from sympy import *
a, b, c, x = symbols("a b c x")
expr = a*x**2 + b*x + c
diff(expr, a), diff(expr, b), diff(expr, c)

7, trong đó

from sympy import *
a, b, c, x = symbols("a b c x")
expr = a*x**2 + b*x + c
diff(expr, a), diff(expr, b), diff(expr, c)

8 là

from sympy import *
a, b, c, x = symbols("a b c x")
expr = a*x**2 + b*x + c
diff(expr, a), diff(expr, b), diff(expr, c)

8 của bạn và

## functions
def exponential_model(x, a, k):
    return a * np.exp(k * x)

def exponential_model_jac(x, a, k):
    return np.array([np.exp(k * x), a * x * np.exp(k * x)]).T  # transform into column

## fit
res, pcov = curve_fit(
    f=exponential_model,
    jac=exponential_model_jac,
    xdata=degradation_df.time,
    ydata=degradation_df.conc,
    p0=[100, -1e-3],
    bounds=(-np.inf, [np.inf, 0]) # k cannot be greater than 0
)

0 là các tham số khác. Có những chữ ký hợp lệ khác, nhưng hãy bắt đầu với chữ ký đơn giản này

## functions
def exponential_model(x, a, k):
    return a * np.exp(k * x)

def exponential_model_jac(x, a, k):
    return np.array([np.exp(k * x), a * x * np.exp(k * x)]).T  # transform into column

## fit
res, pcov = curve_fit(
    f=exponential_model,
    jac=exponential_model_jac,
    xdata=degradation_df.time,
    ydata=degradation_df.conc,
    p0=[100, -1e-3],
    bounds=(-np.inf, [np.inf, 0]) # k cannot be greater than 0
)

1. Không bắt buộc. Hàm trả về Jacobian của hàm của bạn

from sympy import *
a, b, c, x = symbols("a b c x")
expr = a*x**2 + b*x + c
diff(expr, a), diff(expr, b), diff(expr, c)

6. Đôi khi, nó làm cho việc tối ưu hóa trở nên dễ dàng hơn vì trình tối ưu hóa sẽ không phải tính toán gần đúng độ dốc của hàm của bạn bằng số

## functions
def exponential_model(x, a, k):
    return a * np.exp(k * x)

def exponential_model_jac(x, a, k):
    return np.array([np.exp(k * x), a * x * np.exp(k * x)]).T  # transform into column

## fit
res, pcov = curve_fit(
    f=exponential_model,
    jac=exponential_model_jac,
    xdata=degradation_df.time,
    ydata=degradation_df.conc,
    p0=[100, -1e-3],
    bounds=(-np.inf, [np.inf, 0]) # k cannot be greater than 0
)

3. Các giá trị

from sympy import *
a, b, c, x = symbols("a b c x")
expr = a*x**2 + b*x + c
diff(expr, a), diff(expr, b), diff(expr, c)

8

## functions
def exponential_model(x, a, k):
    return a * np.exp(k * x)

def exponential_model_jac(x, a, k):
    return np.array([np.exp(k * x), a * x * np.exp(k * x)]).T  # transform into column

## fit
res, pcov = curve_fit(
    f=exponential_model,
    jac=exponential_model_jac,
    xdata=degradation_df.time,
    ydata=degradation_df.conc,
    p0=[100, -1e-3],
    bounds=(-np.inf, [np.inf, 0]) # k cannot be greater than 0
)

5. Các giá trị

## functions
def exponential_model(x, a, k):
    return a * np.exp(k * x)

def exponential_model_jac(x, a, k):
    return np.array([np.exp(k * x), a * x * np.exp(k * x)]).T  # transform into column

## fit
res, pcov = curve_fit(
    f=exponential_model,
    jac=exponential_model_jac,
    xdata=degradation_df.time,
    ydata=degradation_df.conc,
    p0=[100, -1e-3],
    bounds=(-np.inf, [np.inf, 0]) # k cannot be greater than 0
)

6

## functions
def exponential_model(x, a, k):
    return a * np.exp(k * x)

def exponential_model_jac(x, a, k):
    return np.array([np.exp(k * x), a * x * np.exp(k * x)]).T  # transform into column

## fit
res, pcov = curve_fit(
    f=exponential_model,
    jac=exponential_model_jac,
    xdata=degradation_df.time,
    ydata=degradation_df.conc,
    p0=[100, -1e-3],
    bounds=(-np.inf, [np.inf, 0]) # k cannot be greater than 0
)

7. Tùy chọn (nhưng rất khuyến khích cho các mô hình phi tuyến tính). Dự đoán ban đầu cho các tham số

## functions
def exponential_model(x, a, k):
    return a * np.exp(k * x)

def exponential_model_jac(x, a, k):
    return np.array([np.exp(k * x), a * x * np.exp(k * x)]).T  # transform into column

## fit
res, pcov = curve_fit(
    f=exponential_model,
    jac=exponential_model_jac,
    xdata=degradation_df.time,
    ydata=degradation_df.conc,
    p0=[100, -1e-3],
    bounds=(-np.inf, [np.inf, 0]) # k cannot be greater than 0
)

8. Không bắt buộc. Các ràng buộc đối với các giá trị có thể có của từng tham số (hữu ích để kiểm soát hành vi chung của quá trình tối ưu hóa)

## functions
def exponential_model(x, a, k):
    return a * np.exp(k * x)

def exponential_model_jac(x, a, k):
    return np.array([np.exp(k * x), a * x * np.exp(k * x)]).T  # transform into column

## fit
res, pcov = curve_fit(
    f=exponential_model,
    jac=exponential_model_jac,
    xdata=degradation_df.time,
    ydata=degradation_df.conc,
    p0=[100, -1e-3],
    bounds=(-np.inf, [np.inf, 0]) # k cannot be greater than 0
)

9. Không bắt buộc. Số lần đánh giá chức năng tối đa (về cơ bản, bao nhiêu lần lặp lại)

from sympy import *
a, b, c, x = symbols("a b c x")
expr = a*x**2 + b*x + c
diff(expr, a), diff(expr, b), diff(expr, c)

5 trả về các tham số từ quá trình tối ưu hóa, đây là thứ xác định mô hình

Các mô hình này cũng yêu cầu kiến thức toán học khá tốt (phép tính và hàm cơ bản), để bạn biết mình đang làm gì. Cuối cùng, tôi cũng sẽ sử dụng rộng rãi công cụ tuyệt vời Desmos để vẽ các phương trình của chúng ta. Mã đầy đủ sẽ có sẵn trong repo Github được liên kết ở cuối bài viết

Đoán ban đầu và Jacobian

Đối với các mô hình phi tuyến tính, đôi khi điều quan trọng là chúng tôi cung cấp cho trình tối ưu hóa dự đoán ban đầu tốt (

## functions
def exponential_model(x, a, k):
    return a * np.exp(k * x)

def exponential_model_jac(x, a, k):
    return np.array([np.exp(k * x), a * x * np.exp(k * x)]).T  # transform into column

## fit
res, pcov = curve_fit(
    f=exponential_model,
    jac=exponential_model_jac,
    xdata=degradation_df.time,
    ydata=degradation_df.conc,
    p0=[100, -1e-3],
    bounds=(-np.inf, [np.inf, 0]) # k cannot be greater than 0
)

7) cho các tham số. Hầu hết các bộ đồ tôi sẽ trình bày sẽ không hoạt động nếu không đoán được các thông số của chúng. Trong R, thư viện drc có các quy trình tích hợp giúp người dùng bằng cách chọn cấu hình tham số ban đầu tốt. Tuy nhiên, tôi không tìm thấy thứ gì đó tương tự trong Python, vì vậy chúng tôi sẽ phải dựa vào trực giác toán học và bằng cách xem xét tập dữ liệu của mình. Tuy nhiên, tôi nghĩ rằng những thói quen này có thể được thực hiện dễ dàng, bạn có thể viết một thói quen tìm kiếm dạng lưới chọn một dự đoán tham số ban đầu và sau đó chuyển nó cho

from sympy import *
a, b, c, x = symbols("a b c x")
expr = a*x**2 + b*x + c
diff(expr, a), diff(expr, b), diff(expr, c)

5 thông qua đối số

## functions
def exponential_model(x, a, k):
    return a * np.exp(k * x)

def exponential_model_jac(x, a, k):
    return np.array([np.exp(k * x), a * x * np.exp(k * x)]).T  # transform into column

## fit
res, pcov = curve_fit(
    f=exponential_model,
    jac=exponential_model_jac,
    xdata=degradation_df.time,
    ydata=degradation_df.conc,
    p0=[100, -1e-3],
    bounds=(-np.inf, [np.inf, 0]) # k cannot be greater than 0
)

7 (nhưng điều đó có gì thú vị?)

Đây là một ví dụ về sự không phù hợp, khi bạn không chuyển tham số

## functions
def exponential_model(x, a, k):
    return a * np.exp(k * x)

def exponential_model_jac(x, a, k):
    return np.array([np.exp(k * x), a * x * np.exp(k * x)]).T  # transform into column

## fit
res, pcov = curve_fit(
    f=exponential_model,
    jac=exponential_model_jac,
    xdata=degradation_df.time,
    ydata=degradation_df.conc,
    p0=[100, -1e-3],
    bounds=(-np.inf, [np.inf, 0]) # k cannot be greater than 0
)

7. Mô hình ở đây là Phân rã theo cấp số nhân, được trình bày trong phần tiếp theo và được trang bị trên dữ liệu phân hủy chất phản ứng hóa học. Dữ liệu được thể hiện bằng màu đỏ/hồng

res, pcov = curve_fit(
    f=exponential_model,
    jac=exponential_model_jac, 
    xdata=degradation_df.time,
    ydata=degradation_df.conc,
)

Trong đoạn mã bên dưới, tôi cũng sẽ cung cấp hàm Jacobian cho từng mô hình. Chúng không cần thiết cho sự phù hợp, vì

from sympy import *
a, b, c, x = symbols("a b c x")
expr = a*x**2 + b*x + c
diff(expr, a), diff(expr, b), diff(expr, c)

5 có thể tính gần đúng độ dốc của hàm của bạn bằng số, nhưng tôi sẽ trình bày chúng ở đây để tham khảo (đôi khi bạn có thể cần xem xét các đạo hàm riêng của mô hình của mình, bạn không bao giờ biết . Jacobian ở đây chỉ đơn giản là một vectơ cột với các đạo hàm riêng của hàm w của bạn. r. t. tất cả các thông số của nó. Tôi sẽ sử dụng

## functions
def exponential_decay_asymptote(x, a, b, k):
    return b + (a - b) * np.exp(k * x)


def exponential_decay_asymptote_jac(x, a, b, k):
    return np.array([np.exp(k * x), 1 - np.exp(k * x), x * (a - b) * np.exp(k * x)]).T


offset = 15  # manual concentration offset for illustration purpose
## fit
res, pcov = curve_fit(
    f=exponential_decay_asymptote,
    jac=exponential_decay_asymptote_jac,
    xdata=degradation_df.time,
    ydata=degradation_df.conc + offset,
    p0=[100, 10, -1e-3],  # putting 10 as a starting point for the lower asymptote
    bounds=(-np.inf, [np.inf, np.inf, 0]),  # k cannot be greater than 0
)
a, b, k = res

6 để tự động phân biệt các biểu thức, nhưng bạn có thể sử dụng Wolfram Alpha hoặc bút và giấy cũ tốt

Đây là một ví dụ về cách tính đạo hàm riêng của một hàm bằng cách sử dụng

## functions
def exponential_decay_asymptote(x, a, b, k):
    return b + (a - b) * np.exp(k * x)


def exponential_decay_asymptote_jac(x, a, b, k):
    return np.array([np.exp(k * x), 1 - np.exp(k * x), x * (a - b) * np.exp(k * x)]).T


offset = 15  # manual concentration offset for illustration purpose
## fit
res, pcov = curve_fit(
    f=exponential_decay_asymptote,
    jac=exponential_decay_asymptote_jac,
    xdata=degradation_df.time,
    ydata=degradation_df.conc + offset,
    p0=[100, 10, -1e-3],  # putting 10 as a starting point for the lower asymptote
    bounds=(-np.inf, [np.inf, np.inf, 0]),  # k cannot be greater than 0
)
a, b, k = res

6

from sympy import *
a, b, c, x = symbols("a b c x")
expr = a*x**2 + b*x + c
diff(expr, a), diff(expr, b), diff(expr, c)

$(x^2, x, 1)$

Mô hình lồi/lõm

QUAN TRỌNG. Không phải tất cả các mô hình này đều phi tuyến tính;

Phân rã theo cấp số nhân

Các mô hình phân rã hàm mũ có nhiều ứng dụng khác nhau cho hóa học (phân rã chất), sinh học, kinh tế lượng, v.v.

Tham số hóa chung

$$Y = a\cdot e^{k\cdot X}$$

giải thích thông số

$a$ là giá trị của $Y$ khi $X = 0$
$k$ là mức tăng hoặc giảm tương đối của $Y$ khi tăng một đơn vị $X$

Ở đây chúng tôi sử dụng dữ liệu về phản ứng hóa học phân hủy chất, tôi. e. khi thời gian trôi qua nồng độ của một số chất giảm khi nó xuống cấp. Chúng tôi muốn tìm tốc độ phân hủy chất theo thời gian. Như chúng ta đã thấy trước đó, nếu chúng ta chỉ chuyển hàm cho

from sympy import *
a, b, c, x = symbols("a b c x")
expr = a*x**2 + b*x + c
diff(expr, a), diff(expr, b), diff(expr, c)

5 thì kết quả sẽ không lý tưởng. Vì vậy, chúng tôi cần kiểm tra dữ liệu và mô hình của mình và đưa ra dự đoán ban đầu tốt cho mô hình

Đây là lý do tại sao điều quan trọng là bạn phải biết chính xác những gì bạn đang cố gắng lập mô hình và những thông số đại diện cho điều gì. Phù hợp thì Dự đoán sẽ không hoạt động =)

Từ việc xem xét dữ liệu và tham số hóa của chúng tôi, chúng tôi có thể đoán rằng

giá trị của $a$ phải gần với

## functions
def exponential_decay_asymptote(x, a, b, k):
    return b + (a - b) * np.exp(k * x)


def exponential_decay_asymptote_jac(x, a, b, k):
    return np.array([np.exp(k * x), 1 - np.exp(k * x), x * (a - b) * np.exp(k * x)]).T


offset = 15  # manual concentration offset for illustration purpose
## fit
res, pcov = curve_fit(
    f=exponential_decay_asymptote,
    jac=exponential_decay_asymptote_jac,
    xdata=degradation_df.time,
    ydata=degradation_df.conc + offset,
    p0=[100, 10, -1e-3],  # putting 10 as a starting point for the lower asymptote
    bounds=(-np.inf, [np.inf, np.inf, 0]),  # k cannot be greater than 0
)
a, b, k = res

9

giá trị của

## functions
def negative_exponential_asymptote(x, a, b, k):
    return b - (b - a) * np.exp(-k * x)


def negative_exponential_asymptote_jac(x, a, b, k):
    return np.array(
        [np.exp(-k * x), 1 - np.exp(-k * x), -x * (a - b) * np.exp(-k * x)]
    ).T


## fit
res, pcov = curve_fit(
    f=negative_exponential_asymptote,
    jac=negative_exponential_asymptote_jac,
    xdata=growth_data_df.time,
    ydata=growth_data_df.length,
    p0=[
        0,
        200,
        1e-3,
    ],  # Y can start at 0 when X=0, the asymptote guess is 200 and the rate guess is 1e-3
    bounds=(0, [np.inf, np.inf, 1]),  # adding constraints
)
a, b, k = res

0 phải là số âm, vì dữ liệu cho thấy hành vi phân rã theo cấp số nhân

Hãy nghĩ về

## functions
def negative_exponential_asymptote(x, a, b, k):
    return b - (b - a) * np.exp(-k * x)


def negative_exponential_asymptote_jac(x, a, b, k):
    return np.array(
        [np.exp(-k * x), 1 - np.exp(-k * x), -x * (a - b) * np.exp(-k * x)]
    ).T


## fit
res, pcov = curve_fit(
    f=negative_exponential_asymptote,
    jac=negative_exponential_asymptote_jac,
    xdata=growth_data_df.time,
    ydata=growth_data_df.length,
    p0=[
        0,
        200,
        1e-3,
    ],  # Y can start at 0 when X=0, the asymptote guess is 200 and the rate guess is 1e-3
    bounds=(0, [np.inf, np.inf, 1]),  # adding constraints
)
a, b, k = res

1 với

## functions
def negative_exponential_asymptote(x, a, b, k):
    return b - (b - a) * np.exp(-k * x)


def negative_exponential_asymptote_jac(x, a, b, k):
    return np.array(
        [np.exp(-k * x), 1 - np.exp(-k * x), -x * (a - b) * np.exp(-k * x)]
    ).T


## fit
res, pcov = curve_fit(
    f=negative_exponential_asymptote,
    jac=negative_exponential_asymptote_jac,
    xdata=growth_data_df.time,
    ydata=growth_data_df.length,
    p0=[
        0,
        200,
        1e-3,
    ],  # Y can start at 0 when X=0, the asymptote guess is 200 and the rate guess is 1e-3
    bounds=(0, [np.inf, np.inf, 1]),  # adding constraints
)
a, b, k = res

2. Trong

## functions
def negative_exponential_asymptote(x, a, b, k):
    return b - (b - a) * np.exp(-k * x)


def negative_exponential_asymptote_jac(x, a, b, k):
    return np.array(
        [np.exp(-k * x), 1 - np.exp(-k * x), -x * (a - b) * np.exp(-k * x)]
    ).T


## fit
res, pcov = curve_fit(
    f=negative_exponential_asymptote,
    jac=negative_exponential_asymptote_jac,
    xdata=growth_data_df.time,
    ydata=growth_data_df.length,
    p0=[
        0,
        200,
        1e-3,
    ],  # Y can start at 0 when X=0, the asymptote guess is 200 and the rate guess is 1e-3
    bounds=(0, [np.inf, np.inf, 1]),  # adding constraints
)
a, b, k = res

3, giá trị của

## functions
def exponential_model(x, a, k):
    return a * np.exp(k * x)

def exponential_model_jac(x, a, k):
    return np.array([np.exp(k * x), a * x * np.exp(k * x)]).T  # transform into column

## fit
res, pcov = curve_fit(
    f=exponential_model,
    jac=exponential_model_jac,
    xdata=degradation_df.time,
    ydata=degradation_df.conc,
    p0=[100, -1e-3],
    bounds=(-np.inf, [np.inf, 0]) # k cannot be greater than 0
)

6 sẽ là

## functions
def negative_exponential_asymptote(x, a, b, k):
    return b - (b - a) * np.exp(-k * x)


def negative_exponential_asymptote_jac(x, a, b, k):
    return np.array(
        [np.exp(-k * x), 1 - np.exp(-k * x), -x * (a - b) * np.exp(-k * x)]
    ).T


## fit
res, pcov = curve_fit(
    f=negative_exponential_asymptote,
    jac=negative_exponential_asymptote_jac,
    xdata=growth_data_df.time,
    ydata=growth_data_df.length,
    p0=[
        0,
        200,
        1e-3,
    ],  # Y can start at 0 when X=0, the asymptote guess is 200 and the rate guess is 1e-3
    bounds=(0, [np.inf, np.inf, 1]),  # adding constraints
)
a, b, k = res

5, dẫn đến một giá trị rất nhỏ (trong khi thực tế, nó sẽ dẫn đến giá trị gần bằng 10). Vì vậy,

## functions
def negative_exponential_asymptote(x, a, b, k):
    return b - (b - a) * np.exp(-k * x)


def negative_exponential_asymptote_jac(x, a, b, k):
    return np.array(
        [np.exp(-k * x), 1 - np.exp(-k * x), -x * (a - b) * np.exp(-k * x)]
    ).T


## fit
res, pcov = curve_fit(
    f=negative_exponential_asymptote,
    jac=negative_exponential_asymptote_jac,
    xdata=growth_data_df.time,
    ydata=growth_data_df.length,
    p0=[
        0,
        200,
        1e-3,
    ],  # Y can start at 0 when X=0, the asymptote guess is 200 and the rate guess is 1e-3
    bounds=(0, [np.inf, np.inf, 1]),  # adding constraints
)
a, b, k = res

0 phải nhỏ hơn nhiều, ví dụ như

## functions
def negative_exponential_asymptote(x, a, b, k):
    return b - (b - a) * np.exp(-k * x)


def negative_exponential_asymptote_jac(x, a, b, k):
    return np.array(
        [np.exp(-k * x), 1 - np.exp(-k * x), -x * (a - b) * np.exp(-k * x)]
    ).T


## fit
res, pcov = curve_fit(
    f=negative_exponential_asymptote,
    jac=negative_exponential_asymptote_jac,
    xdata=growth_data_df.time,
    ydata=growth_data_df.length,
    p0=[
        0,
        200,
        1e-3,
    ],  # Y can start at 0 when X=0, the asymptote guess is 200 and the rate guess is 1e-3
    bounds=(0, [np.inf, np.inf, 1]),  # adding constraints
)
a, b, k = res

7

Sau đó, hãy thử khớp với p0=[100, -1e-3] và xem điều gì sẽ xảy ra

## functions
def exponential_model(x, a, k):
    return a * np.exp(k * x)

def exponential_model_jac(x, a, k):
    return np.array([np.exp(k * x), a * x * np.exp(k * x)]).T  # transform into column

## fit
res, pcov = curve_fit(
    f=exponential_model,
    jac=exponential_model_jac,
    xdata=degradation_df.time,
    ydata=degradation_df.conc,
    p0=[100, -1e-3],
    bounds=(-np.inf, [np.inf, 0]) # k cannot be greater than 0
)

Tốt hơn nhiều

Phân rã theo cấp số nhân với tiệm cận thấp hơn

Tham số hóa chung

$$Y = b + (a - b)\cdot e^{k\cdot X}$$

giải thích thông số

$a$ là giá trị của $Y$ khi $X = 0$
$b$ là đường tiệm cận dưới nơi phân rã sẽ ổn định
$k$ là mức tăng hoặc giảm tương đối của $Y$ khi tăng một đơn vị $X$

## functions
def exponential_decay_asymptote(x, a, b, k):
    return b + (a - b) * np.exp(k * x)


def exponential_decay_asymptote_jac(x, a, b, k):
    return np.array([np.exp(k * x), 1 - np.exp(k * x), x * (a - b) * np.exp(k * x)]).T


offset = 15  # manual concentration offset for illustration purpose
## fit
res, pcov = curve_fit(
    f=exponential_decay_asymptote,
    jac=exponential_decay_asymptote_jac,
    xdata=degradation_df.time,
    ydata=degradation_df.conc + offset,
    p0=[100, 10, -1e-3],  # putting 10 as a starting point for the lower asymptote
    bounds=(-np.inf, [np.inf, np.inf, 0]),  # k cannot be greater than 0
)
a, b, k = res

Mô hình tiệm cận (Số mũ âm)

Về cơ bản ngược lại với mô hình phân rã theo cấp số nhân. Nó có thể được sử dụng để mô tả các hiện tượng trong đó $Y$ tăng trưởng có giới hạn khi $X$ tiến đến vô cùng

Ví dụ: nếu $Y$ là chiều dài của một con cá và $X$ đại diện cho tuổi của nó, thì chúng tôi mong đợi sự tăng trưởng (tăng chiều dài của cá) khi nó già đi, nhưng chiều dài của nó sẽ không tăng liên tục, nó có thể sẽ ổn định . Ví dụ về loài cá này dựa trên định luật von Bertalanffy. Ngoài ra, chúng tôi không mong đợi con cá có chiều dài $0$ khi nó được sinh ra, vì vậy tham số hóa chung của loại mô hình này cũng có một tham số để kiểm soát điều này.

Tham số hóa chung

$$Y = b - (b - a)\cdot e^{-k\cdot X}$$

giải thích thông số

$a$ là giá trị của $Y$ khi $X = 0$
$b$ là đường tiệm cận trên nơi tăng trưởng sẽ ổn định (giá trị tối đa có thể đạt được cho $Y$)
$k$ là mức tăng hoặc giảm tương đối của $Y$ khi tăng một đơn vị $X$

## functions
def negative_exponential_asymptote(x, a, b, k):
    return b - (b - a) * np.exp(-k * x)


def negative_exponential_asymptote_jac(x, a, b, k):
    return np.array(
        [np.exp(-k * x), 1 - np.exp(-k * x), -x * (a - b) * np.exp(-k * x)]
    ).T


## fit
res, pcov = curve_fit(
    f=negative_exponential_asymptote,
    jac=negative_exponential_asymptote_jac,
    xdata=growth_data_df.time,
    ydata=growth_data_df.length,
    p0=[
        0,
        200,
        1e-3,
    ],  # Y can start at 0 when X=0, the asymptote guess is 200 and the rate guess is 1e-3
    bounds=(0, [np.inf, np.inf, 1]),  # adding constraints
)
a, b, k = res

Mô hình tiệm cận (bị ràng buộc. bắt đầu từ 0)

Đôi khi mô hình của bạn cần bắt đầu từ 0 (i. e. $a = 0$ từ phương trình trước). Điều này dẫn đến việc tham số hóa đơn giản hơn, mặc dù có cùng hành vi

Tham số hóa chung

$$Y = b\cdot (1 - e^{-k\cdot X})$$

giải thích thông số

$b$ là đường tiệm cận trên nơi tăng trưởng sẽ ổn định (giá trị tối đa có thể đạt được cho $Y$)
$k$ là mức tăng hoặc giảm tương đối của $Y$ khi tăng một đơn vị $X$

## functions
def negative_exponential(x, b, k):
    return b * (1 - np.exp(-k * x))


def negative_exponential_jac(x, b, k):
    return np.array([1 - np.exp(-k * x), b * x * np.exp(-k * x)]).T


offset = -20
## fit
res, pcov = curve_fit(
    f=negative_exponential,
    jac=negative_exponential_jac,
    xdata=growth_data_df.time,
    ydata=growth_data_df.length + offset,
    p0=[200, 1e-3],  # the asymptote guess is 200 and the rate guess is 1e-3
    bounds=(0, [np.inf, 3]),  # adding constraints
)
b, k = res

hồi quy sức mạnh

Tham số hóa chung

$$Y = a\cdot X^{b}$$

giải thích thông số

$a$ là giá trị của $Y$ khi $X = 0$
$b$ là sức mạnh (kiểm soát cách $Y$ tăng so với $X$

Hồi quy lũy thừa tương đương với một đường cong hàm mũ nhưng với logarit của $X$

$$aX^b = a\cdot e^{log(X^b)} = a\cdot e^{b\cdot log(X)}$$

## functions
def power_regr(x, a, b):
    return a * np.power(x, b)


def power_regr_jac(x, a, b):
    return np.array([np.power(x, b), a * np.power(x, b) * np.log(x)]).T


## fit
res, pcov = curve_fit(
    f=power_regr,
    jac=power_regr_jac,
    xdata=species_df.area,
    ydata=species_df.numspecies,
    p0=[200, 2],  # it starts at 0, power guess is 2 (quadratic behavior)
)
a, b = res

Ở đây, tôi đang điều chỉnh mô hình Sức mạnh cho bộ dữ liệu có số lượng loài thực vật theo khu vực lấy mẫu của một số thí nghiệm

khu vực. tổng diện tích của không gian lấy mẫu để đo thực vật
số loài. số loài được tìm thấy

Một thuộc tính hữu ích của tham số Power là bạn có thể ràng buộc nó sao cho phù hợp tùy thuộc vào hành vi tăng trưởng mà bạn đang cố gắng mô hình hóa

$0 < b < 1$. Hình lồi, $Y$ tăng khi $X$ tăng
$b < 0$. Hình lõm, $Y$ giảm khi $X$ tăng
$b > 1$. Hình lõm, $Y$ tăng khi $X$ tăng

Dưới đây là một số ví dụ

Đường cong sigmoidal

Đây là hai đường cong chính để kiểm tra khi bạn cần lập mô hình các sự kiện có hình dạng sigmoidal

đường cong hậu cần

Đường cong logistic phải quen thuộc với bất kỳ nhà khoa học dữ liệu nào. Nó bắt nguồn từ hàm phân phối logistic tích lũy và đối xứng quanh một điểm uốn

Tham số hóa chung

$$Y = a + \frac{c - a}{1 + e^{b(X - d)}}$$

giải thích thông số

$c$ là tiệm cận trên
$a$ là tiệm cận dưới
$d$ kiểm soát vị trí của điểm uốn (so với X)
$b$ độ dốc xung quanh điểm uốn và chúng ta có thể sử dụng nó để kiểm soát hành vi 'tăng trưởng'. Nếu $b < 0$ thì $Y$ tăng khi $X$ tăng và nếu $b > 0$ thì $Y$ giảm khi $X$ tăng

Tham số hóa ở trên là chung, nhưng tùy thuộc vào những gì bạn đang cố gắng điều chỉnh, nó có thể được giảm bớt. Ví dụ: nếu chúng ta thấy $a = 0$, chúng ta có một phương trình đơn giản hơn. Chúng ta cũng có thể buộc $c = 1$ hơn nữa, làm cho phương trình trở nên đơn giản hơn. Tất cả phụ thuộc vào những hạn chế của vấn đề của bạn

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


def logistic_curve_jac(x, a, b, c, d):
    return np.array(
        [
            1 - (1 / (np.exp(b * (-d + x)) + 1)),
            ((-a + c) * (-d + x) * np.exp(b * (-d + x)))
            / (np.exp(b * (-d + x)) + 1) ** 2,
            (1) / (np.exp(b * (-d + x)) + 1),
            (b * (-a + c) * np.exp(b * (-d + x))) / (np.exp(b * (-d + x)) + 1) ** 2,
        ]
    ).T


## fit
res, pcov = curve_fit(
    f=logistic_curve,
    jac=logistic_curve_jac,
    xdata=beetgrowth_df.dae,
    ydata=beetgrowth_df.weightfree,
    p0=[
        0,
        -0.1,
        50,
        30,
    ],  # lower asymptote guess is 0, b is negative (Y increases with X), upper asymptote guess is 50 and inflection point guess is at 30 dae
    bounds=([0, -np.inf, -np.inf, -np.inf], [0.1, 0, np.inf, np.inf]),
    maxfev=3000,
)
a, b, c, d = res

Trong ví dụ này, chúng tôi điều chỉnh Đường cong Logistic để lập mô hình Tăng trưởng Thực vật, trong trường hợp này là Củ cải đường. Các biến là

dae. ngày sau khi cây mọc
không trọng lượng. Đo trọng lượng (gam trên mét vuông chất khô)

Bộ dữ liệu đã được xuất bản trong bài báo này. Scott, R. k. , Wilcockson, S. J. , & Moisey, F. r. (1979). Ảnh hưởng của thời gian trừ cỏ đến sinh trưởng và năng suất củ cải đường

Chức năng Gompertz

Hay cụ thể hơn, Đường cong Gompertz. Không phải lúc nào vấn đề chúng tôi đang lập mô hình cũng đối xứng và hàm Gompertz có thể mô hình hóa các “tốc độ tăng trưởng” khác nhau xung quanh điểm uốn. Tôi khuyên bạn nên xem trang Wikipedia cho chức năng này, đây là một số ví dụ về nơi nó có thể được sử dụng

Mức độ phổ biến của điện thoại di động, trong đó chi phí ban đầu cao (do đó mức độ tiêu thụ chậm), sau đó là giai đoạn tăng trưởng nhanh, tiếp theo là mức tiêu thụ chậm lại khi đạt đến mức bão hòa
Dân số trong một không gian hạn chế, khi tỷ lệ sinh đầu tiên tăng và sau đó chậm lại khi đạt đến giới hạn tài nguyên
Mô hình hóa sự phát triển của khối u
Mô hình hóa tác động thị trường trong tài chính và các khoản vay địa phương tổng hợp năng động

Tham số hóa chung

$$Y = a + (c - a)\cdot e^{-e^{b(X - d)}}$$

Các tham số có cách hiểu tương tự như trong Đường cong Logistic

giải thích thông số

$c$ là tiệm cận trên
$a$ là tiệm cận dưới
$d$ kiểm soát vị trí của điểm uốn
$b$ độ dốc xung quanh điểm uốn và chúng ta có thể sử dụng nó để kiểm soát hành vi 'tăng trưởng'. Nếu $b < 0$ thì $Y$ tăng khi $X$ tăng và nếu $b > 0$ thì $Y$ giảm khi $X$ tăng

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


def gompertz_function_jac(x, a, b, c, d):
    return np.array(
        [
            1 - np.exp(-np.exp(b * (-d + x))),
            -(-a + c) * (-d + x) * np.exp(b * (-d + x)) * np.exp(-np.exp(b * (-d + x))),
            np.exp(-np.exp(b * (-d + x))),
            b * (-a + c) * np.exp(b * (-d + x) * np.exp(-np.exp(b * (-d + x)))),
        ]
    ).T


## fit
res, pcov = curve_fit(
    f=gompertz_function,
    jac=gompertz_function_jac,
    xdata=beetgrowth_df.dae,
    ydata=beetgrowth_df.weightfree,
    p0=[
        0,
        -0.1,
        50,
        30,
    ],  # lower asymptote guess is 0, b is negative (Y increases with X), upper asymptote guess is 50 and inflection point guess is at 30 dae
    bounds=([0, -np.inf, -np.inf, -np.inf], [0.1, 0, np.inf, np.inf]),
    maxfev=3000,
)
a, b, c, d = res

Đường cong Gompertz có một tham số hóa thay thế đảo ngược mô hình tăng trưởng so với mô hình chung

$$Y = a + (c - a)\cdot \left[1 - e^{-e^{-b(X - d)}}\right]$$

Tôi sẽ không phù hợp với bất cứ thứ gì vì nó khá giống với cái trước

Dưới đây là một ví dụ để hình dung rõ hơn về sự khác biệt giữa các đường cong sygmoidal. Các thông số ở đây là.

## functions
def negative_exponential_asymptote(x, a, b, k):
    return b - (b - a) * np.exp(-k * x)


def negative_exponential_asymptote_jac(x, a, b, k):
    return np.array(
        [np.exp(-k * x), 1 - np.exp(-k * x), -x * (a - b) * np.exp(-k * x)]
    ).T


## fit
res, pcov = curve_fit(
    f=negative_exponential_asymptote,
    jac=negative_exponential_asymptote_jac,
    xdata=growth_data_df.time,
    ydata=growth_data_df.length,
    p0=[
        0,
        200,
        1e-3,
    ],  # Y can start at 0 when X=0, the asymptote guess is 200 and the rate guess is 1e-3
    bounds=(0, [np.inf, np.inf, 1]),  # adding constraints
)
a, b, k = res

8

def reverse_gompertz_function(x, a, b, c, d):
    return a + (c - a) * (1 - np.exp(-np.exp(-b * (x - d))))


a, b, c, d = 2, -0.5, 10, 7

Kết luận + Mã

Tôi hy vọng những đoạn mã và ví dụ mã này sẽ hữu ích cho bạn (và cho chính tôi trong tương lai). Có rất nhiều mô hình phi tuyến tính khác mà tôi không đưa vào, một số trong số đó là

phương trình Michaelis-Menten
Đường cong năng suất/mật độ
đường cong log-hậu cần
đường cong Weibull

Như đã đề cập trước đó, bài viết này dựa trên đó bao gồm nhiều chi tiết hơn về các phương trình trên, mặc dù tất cả mã phù hợp đều nằm trong R. Một số bộ dữ liệu tôi đã sử dụng cũng được lấy từ gói

## functions
def negative_exponential_asymptote(x, a, b, k):
    return b - (b - a) * np.exp(-k * x)


def negative_exponential_asymptote_jac(x, a, b, k):
    return np.array(
        [np.exp(-k * x), 1 - np.exp(-k * x), -x * (a - b) * np.exp(-k * x)]
    ).T


## fit
res, pcov = curve_fit(
    f=negative_exponential_asymptote,
    jac=negative_exponential_asymptote_jac,
    xdata=growth_data_df.time,
    ydata=growth_data_df.length,
    p0=[
        0,
        200,
        1e-3,
    ],  # Y can start at 0 when X=0, the asymptote guess is 200 and the rate guess is 1e-3
    bounds=(0, [np.inf, np.inf, 1]),  # adding constraints
)
a, b, k = res

9 trong R

không phải là gì

Hồi quy phi tuyến tính là một loại hồi quy đa thức . Đó là một phương pháp để mô hình hóa mối quan hệ phi tuyến tính giữa các biến phụ thuộc và biến độc lập. Nó được sử dụng tại chỗ khi dữ liệu hiển thị xu hướng cong và hồi quy tuyến tính sẽ không tạo ra kết quả rất chính xác khi so sánh với hồi quy phi tuyến tính.

Python hồi quy tuyến tính và phi tuyến tính là gì?

Một phương trình hồi quy tuyến tính chỉ đơn giản là tính tổng các số hạng. Mặc dù mô hình phải tuyến tính trong các tham số, nhưng bạn có thể tăng một biến độc lập theo số mũ để khớp với một đường cong. Chẳng hạn, bạn có thể bao gồm một thuật ngữ bình phương hoặc lập phương. Các mô hình hồi quy phi tuyến tính là bất kỳ thứ gì không tuân theo một dạng này

không phải là gì

Trong thống kê, hồi quy phi tuyến tính là một dạng phân tích hồi quy trong đó dữ liệu quan sát được mô hình hóa bởi một hàm là sự kết hợp phi tuyến tính của các tham số mô hình và phụ thuộc vào một hoặc nhiều tham số độc lập . Dữ liệu được trang bị bằng một phương pháp xấp xỉ liên tiếp. . The data are fitted by a method of successive approximations.

không phải là gì

plt. show() Hồi quy phi tuyến tính là mối quan hệ giữa các biến độc lập 𝑥 và biến phụ thuộc 𝑦 dẫn đến dữ liệu được mô hình hóa hàm phi tuyến tính . Về cơ bản, bất kỳ mối quan hệ nào không tuyến tính đều có thể được gọi là phi tuyến tính và thường được biểu diễn bằng đa thức 𝑘 độ (lũy thừa cực đại của 𝑥).

programming python

Hồi quy phi tuyến tính Python

Tuyến tính x Phi tuyến

Lắp đường cong trong Python

Đoán ban đầu và Jacobian

Mô hình lồi/lõm

Phân rã theo cấp số nhân

Phân rã theo cấp số nhân với tiệm cận thấp hơn

Mô hình tiệm cận (Số mũ âm)

Mô hình tiệm cận (bị ràng buộc. bắt đầu từ 0)

hồi quy sức mạnh

Đường cong sigmoidal

đường cong hậu cần

Chức năng Gompertz

Kết luận + Mã

không phải là gì

Python hồi quy tuyến tính và phi tuyến tính là gì?

không phải là gì

không phải là gì

Bài Viết Liên Quan

Quảng Cáo

Có thể bạn quan tâm

Toplist được quan tâm

Quảng cáo

Xem Nhiều

Quảng cáo

Chúng tôi

Điều khoản

Trợ giúp

Mạng xã hội