我有一个带有多个峰值的离散数据集 . 我正在尝试生成一种自动方法,用于将高斯曲线拟合到未知数量的数据点 . 最终目标是使用最佳拟合高斯曲线的西格玛值,提供y轴上峰值位置(x轴)的不确定性度量 . 完整的数据集有六个左右的各种形状的独特峰 .

这是一个示例数据集 .

working <- data.frame(age = seq(1, 50), likelihood = c())
likelihood = c(10, 10, 10, 10, 10, 12, 14, 16, 17, 18, 
           19, 20, 19, 18, 17, 16, 14, 12, 11, 10,
           10, 9, 8, 8, 8, 8, 7, 6, 6, 6))

这是高斯拟合程序 . 我在SO上找到了它,但我找不到我再次拿到它的页面,所以请原谅缺乏链接和引用 .

fitG =
function(x,y,mu,sig,scale)
f = function(p){
  d = p[3] * dnorm( x, mean = p[ 1 ], sd = p[ 2 ] )
  sum( ( d - y ) ^ 2)
  }
optim( c( mu, sig, scale ), f )
}

如果我预先定义适合的区域,这很有效 . 例如,仅采用峰值周围的区域并使用输入 mean = 10, sigma = 5和 scale = 1:

work2       <- work[5:20, ]
fit1        <- fitG(work2$age, work2$likelihood, 10, 5, 1)
fitpar1     <- fit1$par
plot(work2$age, work2$likelihood, pch = 20)
lines(work2$age, fitpar1[3]*dnorm(work2$age, fitpar1[1], fitpar1[2]))

但是,我有兴趣以某种方式自动化该过程,我使用 cardidates 包中的 peakwindow 定义整个数据集的峰值中心 . 然后,理想函数将迭代在给定峰值周围拟合中使用的数据点的数量,以便优化高斯参数 . 这是我的尝试:

fitG.2 <- function (x, y) {
  g <- function (z) {
    newdata <- x[(y - 1 - z) : (y + 1 + z), ]
    newfit  <- fitG( newdata$age, newdata$likelihood, 10, 5, 1)
  }
  optimize( f = g, interval = c(seq(1, 100)))
}

但是,我可以't get this type of function to actually work (an error I can'解决) . 我还尝试使用 for 循环创建 function 并设置 break 参数,但此方法对于形状参数变化很大的峰值不一致 . 我可能还有许多其他R功能无法完成此操作 .