R：将数据点稳健地拟合成高斯函数

Question

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

我有一堆（x，y）数据，我想适合Gaussian（又名普通）函数。 重点是，我想要删除的老板。正如人们在下面的示例图中可以看到的那样，还有另一个数据分布会污染我的右边的数据，我不想将它考虑在内来进行拟合（即找到\ sigma，\ mu和整体比例参数）。

[R似乎是这个职位的合适的工具，我发现了一些包（robust，robustbase，MASS例如），它们与强大的配件。

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

我检查了这个related question（和链接的），然而他们把一个单一的向量值作为输入，我有一个向量对，所以我不知道如何转置。

任何帮助如何做到这一点，将不胜感激。

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 

正如指出的那样，对初始值的敏感性总是带着这样的曲线拟合的东西有问题。

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