# 生成一些仿真数据set.seed x <- seqy <- *x)/+x)) + rnorm# 对于一些简单的模型,nls函数可以自动找到合适的参数初值m <- nls)# 计算模型的拟合优度cor)[1] 0.9496598# 将结果可视化plot lines, lty = 2, col = "red", lwd = 3)
输出的图片如下:
library# 利用逻辑斯蒂模型生成人口增长的仿真数据,并用nls估计参数log_growth <- function { with), { dN <- R*N* return)) })}# 逻辑斯蒂增长的参数pars <- c# 设定初值N_ini <- c# 常微分方程的时间阶段times <- seq# 常微分方程out <- ode# 添加一些随机波动N_obs <- out[, 2]+rnorm# 个体数值不能小于1N_obs <- ifelse# 画图plot
这部分代码只是生成了带有随机误差的仿真数据,接下来的部分会展现估计参数初值的技巧。R语言中有一个估计逻辑斯蒂方程参数的内建函数,但它使用的是如下方程:
# 寻找方程的参数SS <- getInitial, data = data.frame)
我们可使用getInitial函数来对模型参数做一个基于数据的初步估计。然后把该函数的输出作为一个向量化参数传递给自启动函数,同时也将无引号的三个参数名赋值给逻辑斯蒂方程。然而,由于SSlogis的参数设定有些不同,我们需要对SSlogis的输出值做一些处理,使得其与逻辑斯蒂方程中的形式一致。
# 改变参数形式K_start <- SS["alpha"]R_start <- 1/SS["scale"]N0_start <- SS["alpha"]/+1)# 构建模型的公式log_formula <- formula/ - 1)))# 拟合模型m <- nls)# 估计参数summaryFormula: N_obs ~ K * N0 * exp/ - 1))Parameters: Estimate Std. Error t value Pr K 1.012e+03 3.446e+01 29.366 <2e-16 ***R 2.010e-01 1.504e-02 13.360 <2e-16 ***N0 9.600e-01 4.582e-01 2.095 0.0415 * ---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 49.01 on 48 degrees of freedomNumber of iterations to convergence: 1 Achieved convergence tolerance: 1.537e-06# 计算拟合优度cor)[1] 0.9910316# 结果可视化lines, col = "red", lty = 2, lwd = 3)
输出图形如下: