R：数据点与高斯函数的稳健拟合答案

【问题标题】：R: Robust fitting of data points to a Gaussian functionR：数据点与高斯函数的稳健拟合
【发布时间】：2013-03-30 17:24:15
【问题描述】：

我需要做一些稳健的数据拟合操作。

我有一堆 (x,y) 数据，我希望它们适合 Gaussian（又名正常）函数。关键是，我想删除 ouliers。正如在下面的示例图中可以看到的那样，还有另一种数据分布污染了右边的数据，我不想考虑它来进行拟合（即找到 \sigma、\mu 和整体尺度参数）。

R 似乎是适合这项工作的工具，我发现了一些与稳健拟合相关的软件包（例如robust、robustbase、MASS）。

但是，他们假设用户已经对 R 有很强的了解，而我的情况并非如此，并且文档仅作为参考手册提供，没有教程或同等内容。我的统计背景相当低，我试图阅读reference material on fitting with R，但它并没有真正帮助（我什至不确定这是正确的方法）。但是我感觉这其实是一个很简单的操作。

我已经检查了这个related question（和链接的），但是它们将单个值向量作为输入，并且我有一个成对向量，所以我不知道如何转置。

如能提供任何帮助，我们将不胜感激。

【问题讨论】：

我认为相关问题是关于将分布拟合为一维数据的密度。你得到的是数据 {x,f(x)} 并且你想拟合 f(x) 的参数，而不是估计分布的参数。
您要删除异常值还是只拟合高斯？
我也有点担心您的数据点看起来不像它们有独立的错误 - 似乎是四个或五个独立的系列。你应该在你的方法中考虑到这一点......
@e4e5f4 ：我想在不考虑异常值的情况下获取底层高斯的参数，所以对我来说这并不重要：要么删除它们，然后计算参数（直接在那个情况下），要么使用一些稳健的拟合算法。
@Spacedman ，第一条评论：是的，嗯......随你喜欢 ;-) 但是知道我该怎么做吗？使用 read.table（2 列）将数据加载到 R 内存中......我被卡住了！

标签： r data-fitting

【解决方案1】：

对数据拟合一条高斯曲线，原理是最小化拟合曲线与数据的平方和差，所以我们定义f我们的目标函数并在其上运行optim：

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)
 }

现在，将其扩展到两个高斯：

fit2G <- function(x,y,mu1,sig1,scale1,mu2,sig2,scale2,...){

  f = function(p){
    d = p[3]*dnorm(x,mean=p[1],sd=p[2]) + p[6]*dnorm(x,mean=p[4],sd=p[5])
    sum((d-y)^2)
  }
  optim(c(mu1,sig1,scale1,mu2,sig2,scale2),f,...)
}

使用第一次拟合的初始参数进行拟合，并对第二个峰值进行目测。需要增加最大迭代次数：

> fit2P = fit2G(data$V3,data$V6,6,.6,.02,8.3,0.10,.002,control=list(maxit=10000))
Warning messages:
1: In dnorm(x, mean = p[1], sd = p[2]) : NaNs produced
2: In dnorm(x, mean = p[4], sd = p[5]) : NaNs produced
3: In dnorm(x, mean = p[4], sd = p[5]) : NaNs produced
> fit2P
$par
[1] 6.035610393 0.653149616 0.023744876 8.317215066 0.107767881 0.002055287

这一切看起来像什么？

> plot(data$V3,data$V6)
> p = fit2P$par
> lines(data$V3,p[3]*dnorm(data$V3,p[1],p[2]))
> lines(data$V3,p[6]*dnorm(data$V3,p[4],p[5]),col=2)

但是我会警惕关于你的函数参数的统计推断......

产生的警告消息可能是由于 sd 参数变为负数。您可以通过使用 L-BFGS-B 并设置下限来解决此问题并获得更快的收敛：

> fit2P = fit2G(data$V3,data$V6,6,.6,.02,8.3,0.10,.002,control=list(maxit=10000),method="L-BFGS-B",lower=c(0,0,0,0,0,0))
> fit2P
$par
[1] 6.03564202 0.65302676 0.02374196 8.31424025 0.11117534 0.00208724

正如所指出的，对初始值的敏感性始终是这样的曲线拟合问题。

【讨论】：

太棒了！这正是我想要的（甚至更多，因为它还给出了“噪声”参数）。我不完全了解所有“R”步骤，但我会详细研究，非常感谢您提供如此清晰准确的答案！我怀疑我会在几周前实现这一目标，也非常感谢。
再强调一点（供以后的读者参考），这种方法对给定的初始值相当敏感，因此必须寻找一个初步的近似值。

【解决方案2】：

拟合高斯：

# your data
set.seed(0)
data <- c(rnorm(100,0,1), 10, 11) 

# find & remove outliers
outliers <- boxplot(data)$out
data <- setdiff(data, outliers)

# fitting a Gaussian
mu <- mean(data)
sigma <- sd(data)

# testing the fit, check the p-value
reference.data <- rnorm(length(data), mu, sigma)
ks.test(reference.data, data)

【讨论】：