R语言实践之三次函数批量拟合及积分计算

最新推荐文章于 2023-08-18 11:14:07 发布

阿丘的地理空间大喇叭

最新推荐文章于 2023-08-18 11:14:07 发布

阅读量1.7k

点赞数 1

文章标签： r语言开发语言

本文链接：https://blog.csdn.net/weixin_44122282/article/details/124758687

版权

R语言实践之三次函数批量拟合及积分计算

前段时间写论文时数据处理阶段遇到一个问题，需要拟合三次函数并计算峰值及对应的积分面积(如下图)，而且数据量比较大，今天选取部分数据展示一下操作过程。
图1. 三次函数拟合曲线图解
先看一下数据，一共有15个sheet表，每个sheet里面包括ID、0-1000的样本点，其中还有一些NA值需要剔除，总共1527组数据需要拟合。
图2.数据截图
首先读入数据，用filter把小于0的数据剔除掉。先拟合一组数据试试，提取第一行命名为y，再添加一列x值，把格式转为数据框(为后面调用lm做准备)

  f<-"adata.xlsx"
  data <- read_excel(f,sheet=1)
  #去除空值数据
  data_lm <- filter(data,`0`>0 ,`100`>0,`200`>0,`300`>0,`400`>0,
                    `500`>0,`600`>0,`700`>0,`800`>0,`900`>0,`1000`>0)
  data_lmd <- t(select(data_lm[1,],2:12))#提取第一行
  datan <- as.data.frame(cbind(x,data_lmd))
  colnames(datan) <- c("x","y")

再来看看怎么拟合，那第一步就建立一个function，里面包含要拟合的三次函数的模型，公式就是y~a+bx+cx2+d*x3

lm3c <- function(x,y,datan){
fit2<-lm( y ~ x+ I(x^2)+I(x^3), data = datan )
  myfun <- function(x){
    p <-fit2$coefficients[1]+fit2$coefficients[2]*x+fit2$coefficients[3]*x^2+fit2$coefficients[4]*x^3
    return (p)
  }
}

三次函数已经建立好了，那求其峰值也就是求一阶导函数的零点，所以建立其导函数myfun_dao

###导函数
  myfun_dao <- function(y){
    f <- fit2$coefficients[2]+2*fit2$coefficients[3]*y+3*fit2$coefficients[4]*y^2
    return (f)
  }

在这过程中又发现，有的数据拟合出来其一阶导没有零点，但实际情况确实有逐渐放缓的趋势(图3)。所以就建立二阶导myfun_dao2，如果一阶导没有零点，就判断二阶导，其与x轴的交点即为放缓趋势的临界点。

  ###二阶导数函数
  myfun_dao2 <- function(y){
    f <- 2*fit2$coefficients[3]+6*fit2$coefficients[4]*y
    return (f)
  }

在这里插入图片描述

在这里插入图片描述
图3.一阶导无零点情况(a:原函数; b:导函数)
接下来到了计算导函数的根，R里面的函数为uniroot，它从给定的(x,y)区间查找函数f的根，这个过程比较笨，前提是需要你给出的x到y的区间内一定要包含根，不然算不出来。这就意味着你要先确定根的范围，那这不进入了死循环？我都知道根的范围了我还算啥呢。通过help查询uniroot的解释说明，提到可以设定extendInt参数，来扩展它的端点值，但通过验证，实际上这个扩展是有限的，只能是扩展一定范围，而不是整个x轴。所以在我的数据里面x范围为0-1000，，如下面的代码所示，我以100为间隔划分了10个区间进行迭代寻根，暂时没有发现报错。

在这里插入图片描述

#划定范围寻根
  troot <- try(uniroot(myfun_dao,c(0,100), extendInt ="yes",tol=0.0001)$root)
  for(k in seq(200,1000,by=100)){
    #%in%判断前面一个向量内的元素是否在后面一个向量中，返回布尔值
    if("try-error" %in% class(try(uniroot(myfun_dao,c(0,k), extendInt ="yes",tol=0.0001)$root))){next} 
    else{textent <- try(uniroot(myfun_dao,c(0,k), extendInt ="yes",tol=0.0001)$root)}
  }
  if(textent==1000){
    textent <- uniroot(myfun_dao2,c(0,1000), extendInt ="yes",tol=0.0001)$root
  }

找到根也就是峰值之后，后面就好办了。用integrate函数计算积分面积，最后返回峰值和积分面积两个参数。

    extent1 <- textent
    extent_integ <- integrate(myfun,0,extent1)$value##算积分面积
    ep <- c(extent1,extent_integ)##组合返回变量
    return (ep)

把这一套拟合的方法封装到一个工具里面，以便调用，完整的function如下，命名lm3c

lm3c <- function(x,y,datan){
  ##定义拟合函数
  fit2<-lm( y ~ x+ I(x^2)+I(x^3), data = datan )
  myfun <- function(x){
    p <- fit2$coefficients[1]+fit2$coefficients[2]*x+fit2$coefficients[3]*x^2+fit2$coefficients[4]*x^3
    return (p)
  }
  ###导数函数
  myfun_dao <- function(y){
    f <- fit2$coefficients[2]+2*fit2$coefficients[3]*y+3*fit2$coefficients[4]*y^2
    return (f)
  }
  ###二阶导数函数
  myfun_dao2 <- function(y){
    f <- 2*fit2$coefficients[3]+6*fit2$coefficients[4]*y
    return (f)
  }
  #划定范围寻根
  troot <- try(uniroot(myfun_dao,c(0,100), extendInt ="yes",tol=0.0001)$root)
  for(k in seq(200,1000,by=100)){
    #%in%判断前面一个向量内的元素是否在后面一个向量中，返回布尔值
    if("try-error" %in% class(try(uniroot(myfun_dao,c(0,k), extendInt ="yes",tol=0.0001)$root))){next}
    else{textent <- try(uniroot(myfun_dao,c(0,k), extendInt ="yes",tol=0.0001)$root)}
  }
  if(textent==1000){
    textent <- uniroot(myfun_dao2,c(0,1000), extendInt ="yes",tol=0.0001)$root
  }
    extent1 <- textent
    extent_integ <- integrate(myfun,0,extent1)$value##算积分面积
    ep <- c(extent1,extent_integ)##组合返回变量
    return (ep)
}

最后加入循环迭代，即可拟合出所有的数据，拟合的结果放到list里面

#初始化参数
x <- c(0,100,200,300,400,500,600,700,800,900,1000)
lm11 <- list()
class1_epf <- list()
#读取数据
for (i in 1:15){
  f<-"adata.xlsx"
  data <- read_excel(f,sheet=i)
  #去除空值数据
  data_lm <- filter(data,`0`>0 ,`100`>0,`200`>0,`300`>0,`400`>0,
                    `500`>0,`600`>0,`700`>0,`800`>0,`900`>0,`1000`>0)
  #先拟合第一行，再累加起来
  data_lmd <- t(select(data_lm[1,],2:12))#提取第一行
  datan <- as.data.frame(cbind(x,data_lmd))
  colnames(datan) <- c("x","y")
  lm11 <- lm3c(x,y,datan)#拟合第一行
  lm11 <- as.data.frame(cbind(lm11[[1]],lm11[[2]]))
  names(lm11) <- c("extent","PCI")
  #拟合剩下的行
  for(j in 2:length(data_lm[[1]])){
    lm1 <- list()
    data_lmd <- t(select(data_lm[j,],2:12))#提取单行
    datan <- as.data.frame(cbind(x,data_lmd))
    colnames(datan) <- c("x","y")
    lm1 <- lm3c(x,y,datan)#拟合
    lm1 <- as.data.frame(cbind(lm1[[1]],lm1[[2]]))
    names(lm1) <- c("extent","PCI")
    lm11 <- rbind(lm11,lm1)#合并
  }
  class1_ep <- cbind(select(data_lm,1),lm11)#加上名称、ID 
}