R语言采用优化方法拟合曲线并计算AIC,BIC,LRT

文章目录前言一、RMSD(均方根误差)优化1.RMSD(均方根误差)优化代码2.运行结果二、梯度下降(Gradient Descent)1.梯度下降(Gradient Descent)优化代码2.运行结果前言采用冥函数为模型：Y=x^(-1)采用优化标准和方法：RMSD(均方根误差)，梯度下降(Gradient Descent)对比BIC和AIC一、RMSD(均方根误差)优化1.RMSD(均方根误

椰麦吸

11713人浏览 · 2021-11-28 17:00:52

椰麦吸 · 2021-11-28 17:00:52 发布

文章目录

前言
一、R代码实现
二、运行结果

前言

用R语言(或其他)编写和调试一个拟合记忆任务数据的程序，基于所设计的参数，分别用冥函数和指数函数作为模型，利用多种方法（如RMSD、最大似然、贝叶斯等）进行参数寻优，利用LRT指标来对比相同模型不同自由参数模型的拟合优劣，最后利用多种模型间对比指标（如AIC，BIC等）对两种模型的优劣进行对比。目的是学会运用机器学习的方法来进行简单的数据建模和参数寻优，熟悉简单机器学习算法的开发过程并理解其实现原理

采用指数函数为模型：
优化的目标函数为： Y=2^(-x)

采用优化标准和方法：
RMSD(均方根误差)；
梯度下降(Gradient Descent)；
最大似然估计。

对比BIC和AIC

一、R代码实现

1.导入库

#导入包maxLik，最大使然估计使用
library(maxLik)

注：本文使用的软件是RStudio，包的安装请参考网上教程

2.随机生成原始数据

#随机生成原始数据
sample_count <- 25 # 设定样本容量
x <- runif(sample_count, min = -4, max = 4)
y <- 2^(-x)+rnorm(sample_count,mean=0,sd =0.3)

3.RMSD

#RMSD:
RMSD <- function(parms,y,x){
  if(any(parms<1)|| any(parms>3.5)) return(1e6)
  pow_pred <- (parms["a"])^(-x)
  return(sqrt(sum((pow_pred-y)^2)/length(x)))
}
sparms <- c(3)
names(sparms) <- c("a")
pout <- optim(sparms,RMSD,y=y,x=x)
pow_pred <- (pout$par["a"])^(-x)

4.梯度下降

#梯度下降:
gd <- function(x,y,a1,alpha,tol,M){
  i <- 1
  res <- data.frame(a1=a1)
  repeat{
    J0 <- 1/2*mean((a1^(-x)-y)^2)
    a1_hat <- a1 - alpha*mean((a1^(-x)-y)*(-x)*a1^(-x-1)) #mean()式子是J0对a1求导
    a1 <- a1_hat
    J1 <- 1/2*mean((a1_hat^(-x)-y)^2)
    res <- rbind(res,data.frame(a1=a1))
    if( abs(J1-J0) < tol | i >= M ) 
      break
    i <- i + 1
  }
  return(res)
}
res.gd <- gd(x = x,y = y,a1 = 3,alpha=0.01,tol=1e-8,M=10000)
a_2<- tail(res.gd,1)
a2 <- unlist(rep(a_2,length(x)))
y_pre <- a2^(-x)

5.最大似然估计

#最大使然估计:
logL <- function(params_vector){
  a <- params_vector[1]
  sigma <- params_vector[2]
  sum(dnorm(y-a^(-x), sd = sigma, log = TRUE))
}
start_params <- c(3,1)
result <- maxLik(logLik = logL, start = start_params)
Estimate_a <-summary(result)$estimate[1,1] #提取参数a
logL_pred <- Estimate_a^(-x)

6.做出优化后图像

x11()
par(cex.axis=1.2,cex.lab=1.4)
par(mar=(c(5,5,3,2)+0.1),las=1)
plot(x,y,
     xlab="X",
     ylab="Y",
     ylim=c(0,16),xlim=c(-4,4),xaxt="n",type="n",main = "指数函数优化拟合")
lines(spline(x,pow_pred),lwd=2,col="red")  #RMSD优化拟合曲线
lines(spline(x,y_pre),lwd=2,col="blue")  #梯度下降优化拟合曲线
lines(spline(x,logL_pred),lwd=2,col="green")  #最大使然优化拟合曲线
points(x,y,pch=21,bg="dark grey",cex=2) #画出设定点
axis(1,at=c(-4:4)) #x轴

7.求AIC,BIC

k_1 <- 1
n_1 <-c(length(y))
#根据AIC=2k+nIn(SSR/n) K:参数个数   n: 观察数  SSR:残差平方和
AIC_M <- function(y_test,y_pred,k,n){
  resid <- y_test-y_pred
  SSR <- sum(resid ** 2)
  return(2*k+n*log(SSR/n)) 
}
#根据BIC=KIn(n)+nIn(SSR/n) K:参数个数   n: 观察数  SSR:残差平方和
BIC_M <- function(y_test,y_pred,k,n){
  resid <- y_test-y_pred
  SSR <- sum(resid ** 2)
  return((k*log(n)+n*log(SSR/n))
}

#(RMSD)指数函数AIC
AIC_M(y,pow_pred,k_1,n_1)
#(RMSD)指数函数BIC
BIC_M(y,pow_pred,k_1,n_1)

#(梯度下降)指数函数AIC
AIC_M(y,y_pre,k_1,n_1)
#(梯度下降)指数函数BIC
BIC_M(y,y_pre,k_1,n_1)

#(最大使然)指数函数AIC
AIC_M(y,logL_pred,k_1,n_1)
#(最大使然)指数函数BIC
BIC_M(y,logL_pred,k_1,n_1)

#(RMSD)指数函数优化后参数
pout$par["a"]

#(梯度下降)指数函数优化后参数
a_2[1,1]

#(最大似然估计)指数函数优化后参数
Estimate_a

8.求LRT

h0 <- 2
#LRT评价指标：
LRT <- function(h0,h){
  LRT_pre <- prod(h^(-x))
  LRT_rel <- prod(h0^(-x))
  return(LRT_pre/LRT_rel)
}

#RMSD的LRT：
LRT(h0,pout$par["a"])

#梯度下降的LRT：
LRT(h0,a_2[1,1])

#最大似然的LRT：
LRT(h0,Estimate_a)

二、运行结果

1.图像输出

在这里插入图片描述

由于原始数据是随机生成，会导致图像会有差异，建议多运行几次，找到自己心仪的图像，再者图像原本是三条预测曲线组成，画出后的图像略有重合。最后本人才学疏浅，代码有几处bug，希望你们见谅，但大体能得到所要结果！！！
RMSD优化拟合曲线(红色)
梯度下降优化拟合曲线(蓝色)
最大似然优化拟合曲线(绿色)