R：生成混合分布的函数答案

【问题标题】：R : function to generate a mixture distributionR：生成混合分布的函数
【发布时间】：2014-06-22 05:31:46
【问题描述】：

我需要从混合分布中生成样本

40% 的样本来自 Gaussian(mean=2,sd=8)
20% 的样本来自 Cauchy(location=25,scale=2)
40% 的样本来自 Gaussian(mean = 10, sd=6)

为此，我编写了以下函数：

dmix <- function(x){
prob <- (0.4 * dnorm(x,mean=2,sd=8)) + (0.2 * dcauchy(x,location=25,scale=2)) + (0.4 * dnorm(x,mean=10,sd=6))
return (prob)
}

然后测试：

foo = seq(-5,5,by = 0.01)
vector = NULL
for (i in 1:1000){
vector[i] <- dmix(foo[i])
}
hist(vector)

我得到这样的直方图（我知道这是错误的）-

我做错了什么？谁能指点一下？

【问题讨论】：

我只注意到一件事：您可以使用 hist(dmix(seq(-5,5,by = 0.01))) without looping 和 vector 创建相同的情节
你能做一个随机样本的直方图吗？ dmix <- function(x = 100) {prob <- c(rnorm(x * .4, mean = 2, sd = 8), rcauchy(x * .2, location = 25, scale = 2), rnorm(x * .4, mean = 10, sd = 6)); hist(prob)}; dmix()
@rawr 我可以使用这个 - dmix2

标签： r probability markov-chains mcmc mixture-model

【解决方案1】：

如果可以，请始终使用 R 向量化。即使实际上丢弃了许多值，它通常也更有效。（至少比以前的解决方案更快，并且避免了额外的库）

rmy_ve = function(n){

##generation of (n x 3) matrix. 
##Each column is a random sample of size n from a single component of the mixture
temp = cbind(rnorm(n,2,8),rcauchy(n,25,2),rnorm(n,10,6))

##random generation of the indices
id = sample(1:3,n,rep = T,prob = c(.4,.2,.4))  
id = cbind(1:n,id)
temp[id]
}


> microbenchmark(rmy_ve(1e6),rmyMix(1e6))
Unit: milliseconds
       expr     min       lq     mean   median       uq      max    neval
rmy_ve(1e+06) 241.904 248.7528 272.9119 260.8752 298.1126 322.7429   100
rmyMix(1e+06) 270.917 322.3627 341.4779 329.1706 364.3947 561.2608   100

【讨论】：

【解决方案2】：

当然还有其他方法可以做到这一点，但是 distr 包使它非常简单。（See also this answer 是另一个示例以及关于 distr 和朋友的更多详细信息）。

library(distr)

## Construct the distribution object.
myMix <- UnivarMixingDistribution(Norm(mean=2, sd=8), 
                                  Cauchy(location=25, scale=2),
                                  Norm(mean=10, sd=6),
                                  mixCoeff=c(0.4, 0.2, 0.4))
## ... and then a function for sampling random variates from it
rmyMix <- r(myMix)

## Sample a million random variates, and plot (part of) their histogram
x <- rmyMix(1e6)
hist(x[x>-100 & x<100], breaks=100, col="grey", main="")

如果您只想直接查看混合分布的 pdf，请执行以下操作：

plot(myMix, to.draw.arg="d")

【讨论】：

这似乎是我要找的。span>
哦，这正是我要找的！非常感谢！不知道“发行版”的存在。确实遇到了“mixtools”。然而，这似乎更多的是用于分析而不是生成。
@Raaj 太好了。很高兴这有帮助。正如我确定您刚刚看到的那样，这只是触及了 Ruckdeschel 等人的内容的表面。已经拼好了！
我可以再问一个问题吗？在我看来， [-100 to 100] 类似于截断（或者相同？）它的目的是什么？我当然会查看文档，但如果可以的话，请告诉我。谢谢。
当你绘制柯西随机变量时，如果你不以这种方式剪掉尾巴，那么你最终会得到一个 x 范围非常大的直方图（例如 -1e-7 到 +1e7），然后你根本看不到分布的主体。