2012-11-21 134 views
2

我在R中搜索一个函数/程序包名称,它允许分离两个叠加的正态分布。分布看起来是这样的:在R中分离两个叠加的正态分布

x<-c(3.95, 3.99, 4.0, 4.04, 4.1, 10.9, 11.5, 11.9, 11.7, 12.3) 
+2

为什么这种被关闭作为题外话? – nico

+0

鉴于良好的答案,我加入我的投票重新开放,但我支持我的原始投票。没有足够的细节来说明他们想要实现什么方法来使其成为一个单纯的编程问题,而且它的质量太差,无法迁移到CV。 –

+0

@ AriB.Friedman:...并且根据当前的简历政策,这个问题会被认为与CV相关的程序设计过于相关。 – russellpierce

回答

9

我在过去使用向量广义线性模型有很好的结果。 VGAM package对此很有用。

mix2normal1函数允许估计两个单变量正态分布混合的参数。

小例如

require(VGAM) 
set.seed(12345) 

# Create a binormal distribution with means 10 and 20 
data <- c(rnorm(100, 10, 1.5), rnorm(200, 20, 3)) 

# Initial parameters for minimization algorithm 
# You may want to create some logic to estimate this a priori... not always easy but possible 
# m, m2: Means - s, s2: SDs - w: relative weight of the first distribution (the second is 1-w) 
init.params <- list(m=5, m2=8, s=1, s2=1, w=0.5) 

fit <<- vglm(data ~ 1, mix2normal1(equalsd=FALSE), 
       iphi=init.params$w, imu=init.params$m, imu2=init.params$m2, 
       isd1=init.params$s, isd2=init.params$s2) 

# Calculated parameters 
pars = as.vector(coef(fit)) 
w = logit(pars[1], inverse=TRUE) 
m1 = pars[2] 
sd1 = exp(pars[3]) 
m2 = pars[4] 
sd2 = exp(pars[5]) 

# Plot an histogram of the data 
hist(data, 30, col="black", freq=F) 
# Superimpose the fitted distribution 
x <- seq(0, 30, 0.1) 
points(x, w*dnorm(x, m1, sd1)+(1-w)*dnorm(x,m2,sd2), "l", col="red", lwd=2) 

这正确地给出( “真” 参数 - 10,20,1.5,3)

> m1 
[1] 10.49236 
> m2 
[1] 20.06296 
> sd1 
[1] 1.792519 
> sd2 
[1] 2.877999 

Fit bimodal distribution

3

您可能需要使用nls,非线性回归工具(或其他NONLIN回归量)。我猜你有一个代表叠加分布的数据向量。然后,大致上,nls(y~I(a*exp(-(x-meana)^2/siga) + b*exp(-(x-meanb)^2/sigb)),{initial guess values required for all constants}),其中y是您的分布,x是域。 我根本没有考虑这个问题,所以我不确定哪种收敛方法不太可能失败。