R语言里面，用nls()非线性最小二乘，就在自带的stats包里面。下面是一个例子——
m.lrp <- nls(y ~ f.lrp(x, a, b, t.x), data = test, start = list(a = 0, b = 0.1,
t.x = 50), trace = T, control = list(warnOnly = T, minFactor = 1/2048))
## 32.74 : 0.0 0.1 50.0
## 7.352 : 2.16251 0.04619 59.34899
## 1.25 : 2.16251 0.04619 70.24081
## 1.116 : 2.15689 0.04639 72.09071
## 1.116 : 2.15689 0.04639 72.08250
summary(m.lrp)
##
## Formula: y ~ f.lrp(x, a, b, t.x)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 2.15689 0.06562 32.9 <2e-16 ***
## b 0.04639 0.00157 29.6 <2e-16 ***
## t.x 72.08250 1.76996 40.7 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.176 on 36 degrees of freedom
##
## Number of iterations to convergence: 4
## Achieved convergence tolerance: 3.63e-09
画图来看看拟合的可靠性：
plot(test$y ~ test$x, main = "Linear response and plateau yield response", xlab = "Fertilizer added",
ylab = "Crop yield")
(max.yield <- a.0 + b.0 * t.x.0)
## [1] 5.5
lines(x = c(0, t.x.0, 120), y = c(a.0, max.yield, max.yield), lty = 2, col = "blue")
abline(v = t.x.0, lty = 3, col = "blue")
abline(h = max.yield, lty = 3, col = "blue")
(max.yield <- coefficients(m.lrp)["a"] + coefficients(m.lrp)["b"] * coefficients(m.lrp)["t.x"])
## a
## 5.501
lines(x = c(0, coefficients(m.lrp)["t.x"], 120), y = c(coefficients(m.lrp)["a"],
max.yield, max.yield), lty = 1)
abline(v = coefficients(m.lrp)["t.x"], lty = 4)
abline(h = max.yield, lty = 4)
text(120, 4, "known true model", col = "blue", pos = 2)
text(120, 3.5, "fitted model", col = "black", pos = 2)
Python的话，可以用scipy扩展库的leastsq函数。下面是一个例子——
import numpy as np
from scipy.optimize import leastsq
def residuals(p, x, y):
temp = np.array([0,0,0,0,0,0],dtype=float)
for i in range(0, len(x)):
temp = fun(x, p)
return temp - y
r = leastsq(residuals, [1, 1], args=(x, y))
扫码加好友，拉您进群
全部版块
我的主页
收藏
