以下代码是在R语言中利用nls函数对胸径生长量进行拟合的代码,先不想采用随机抽样法划分建模与检验数据,而是使用十折交叉验证去测试拟合的精度,并获得拟合模型的AIC/BIC和十折交叉的RMSEcv、R2cv、和MAEcv值作为评价指标,请写出完整代码。data <- read.csv("E:/data/Dbaihua.csv",header = T)
set.seed(123) # 设置种子值,保证随机抽样可重复性
rows <- nrow(data) # 计算数据的行数
ns <- sample(1:rows,round(0.75*rows)) # ns为抽出的作为建模数据的行号
ns <- sort(ns) # 将ns按从小到大排序
model.data <- data[ns,] # 建模数据,占总数据的75%
test.data <- data[-ns,]
#### 拟合 t-D ####
# Schumacher 方程
n1 <- nls(D~A*exp(-b/t),
data = model.data,
start = c(A=33,b=24))
#Mitscherlich
n2 <- nls(D~a*(1-b*exp(-c*t)),
data = model.data,
start = c(a=40,b=0.1,c=0.05))
#Logistic方程
n3 <- nls(D~A/(1+m*exp(-r*t)),
data = model.data,
start = c(A=20,m=20,r=0.1))
#Gompertz 方程
n4 <- nls(D~A*exp((-b)*exp(-r*t)),
data = model.data,
start = c(A=20,b=5,r=0.05))
# 考尔夫
n5 <- nls(D~a*exp(-b*t^(-c)),
data = model.data,
start = c(a=108,b=11,c=0.5))
#理查德(Richards)方程
n6 <- nls(D~A*(1-exp(-r*t))^c,
data = model.data,
start = c(A=25,r=0.05,c=1))
# 计算拟合指标 决定系数(R2)
1-sum(resid(n1)^2)/sum((model.data[,"D"]-mean(model.data[,"D"]))^2)
1-sum(resid(n2)^2)/sum((model.data[,"D"]-mean(model.data[,"D"]))^2)
1-sum(resid(n3)^2)/sum((model.data[,"D"]-mean(model.data[,"D"]))^2)
1-sum(resid(n4)^2)/sum((model.data[,"D"]-mean(model.data[,"D"]))^2)
1-sum(resid(n5)^2)/sum((model.data[,"D"]-mean(model.data[,"D"]))^2)
1-sum(resid(n6)^2)/sum((model.data[,"D"]-mean(model.data[,"D"]))^2)
# 计算拟合指标均方根误差(RMSE)
sqrt(sum(resid(n1)^2)/(length(model.data[,"D"])-length(coef(n1))))
sqrt(sum(resid(n2)^2)/(length(model.data[,"D"])-length(coef(n2))))
sqrt(sum(resid(n3)^2)/(length(model.data[,"D"])-length(coef(n3))))
sqrt(sum(resid(n4)^2)/(length(model.data[,"D"])-length(coef(n4))))
sqrt(sum(resid(n5)^2)/(length(model.data[,"D"])-length(coef(n5))))
sqrt(sum(resid(n6)^2)/(length(model.data[,"D"])-length(coef(n6))))
# 计算拟合指标均方误差(MSE)
(sum(resid(n1)^2)/(length(model.data[,"D"])-length(coef(n1))))
(sum(resid(n2)^2)/(length(model.data[,"D"])-length(coef(n2))))
(sum(resid(n3)^2)/(length(model.data[,"D"])-length(coef(n3))))
(sum(resid(n4)^2)/(length(model.data[,"D"])-length(coef(n4))))
(sum(resid(n5)^2)/(length(model.data[,"D"])-length(coef(n5))))
(sum(resid(n6)^2)/(length(model.data[,"D"])-length(coef(n6))))
# 计算AIC
AIC(n1,n2,n3,n4,n5,n6)
BIC(n1,n2,n3,n4,n5,n6)
p1=predict(n1,test.data)
e1=test.data$D-p1
(MARE = mean(abs(e1/test.data$D)))
(MAE = mean(abs(e1)))
p2=predict(n2,test.data)
e2=test.data$D-p2
(MARE = mean(abs(e2/test.data$D)))
(MAE = mean(abs(e2)))
p3=predict(n3,test.data)
e3=test.data$D-p3
(MARE = mean(abs(e3/test.data$D)))
(MAE = mean(abs(e3)))
p4=predict(n4,test.data)
e4=test.data$D-p4
(MARE = mean(abs(e4/test.data$D)))
(MAE = mean(abs(e4)))
p5=predict(n5,test.data)
e5=test.data$D-p5
(MARE = mean(abs(e5/test.data$D)))
(MAE = mean(abs(e5)))
p6=predict(n6,test.data)
e6=test.data$D-p6
(MARE = mean(abs(e6/test.data$D)))
(MAE = mean(abs(e6)))
summary(n1)
summary(n2)
summary(n3)
summary(n4)
summary(n5)
summary(n6)
最新发布