R:将数据点稳健拟合到高斯函数

Question

我需要做一些强大的数据拟合操作.

我有一堆(x,y)数据,我想要适合高斯(又名正常)函数.关键是,我想删除ouliers.正如可以在下面的示例图中看到的那样,还有另一个数据分布在右边污染我的数据,我不想考虑它来进行拟合(即找到\ sigma,\ mu和整体规模参数).

R似乎是适合这项工作的工具,我发现了一些与稳健拟合相关的软件包(健壮,强大的基础,例如MASS).

但是,他们假设用户已经非常了解R,这不是我的情况,文档仅作为一种参考手册提供,没有教程或同等学历.我的统计背景相当低,我试图阅读关于拟合R的参考资料,但它并没有真正帮助(我甚至不确定这是正确的方法).但我觉得这实际上是一个非常简单的操作.

我已经检查了这个相关的问题(以及链接的问题),但是它们将单个矢量值作为输入,并且我有一对矢量,所以我看不到如何转置.

任何有关如何做到这一点的帮助将不胜感激.

Answer 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

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