在R中拟合多模态分布；从拟合的分布生成新值

Question

我正在处理小样本数据：

>dput(dat.demand2050.unique)  
c(79, 56, 69, 61, 53, 73, 72, 86, 75, 68, 74.2, 80, 65.6, 60, 54)    

所对应的密度分布如下图所示：

我知道这些值来自于两个不同的区间 - 低和高，并且假设底层过程服从正态分布，我使用了 mixtools 包来拟合一个双峰分布：

set.seed(99)  
dat.demand2050.mixmdl <- normalmixEM(dat.demand2050.unique, lambda=c(0.3,0.7), mu=c(60,70), k=2)

这使我得到了以下结果：

（实线为拟合曲线，虚线为原始密度）。

# get the parameters of the mixture
dat.demand2050.mixmdl.prop <- dat.demand2050.mixmdl$lambda    #mix proportions
dat.demand2050.mixmdl.means <- dat.demand2050.mixmdl$mu    #modal means
dat.demand2050.mixmdl.dev <- dat.demand2050.mixmdl$sigma   #modal std dev  

混合参数为：

>dat.demand2050.mixmdl.prop  #mix proportions  
[1] 0.2783939 0.7216061  
>dat.demand2050.mixmdl.means  #modal means  
[1] 56.21150 73.08389  
>dat.demand2050.mixmdl.dev  #modal std dev  
[1] 3.098292 6.413906 

我有以下问题：

1. 为了生成一个接近真实分布的新值集合，我的方法正确吗？还是有更好的流程？
2. 如果我的方法正确，如何使用这个结果来生成一组符合该混合分布的随机值？

- avg

---

我认为这个问题可能更适合于CrossValidated：http://stats.stackexchange.com - David Marx

@DavidMarx 是的，我曾经考虑过这个问题，甚至是否要跨帖，但最终决定在这里写，因为我的第二个问题更多涉及编码。然而，如果版主认为那里更合适，我很乐意这样做。 - avg

我不确定你的方法是否明智。你没有说明你打算如何使用随机数。此外，你的样本量非常小，从这么小的样本量估计正态分布有些可疑。也许自助法对于你的最终目标会是更好的方法？ - Roland

@ Roland 确实，样本大小很小，但这就是我所拥有的。数据来自一组研究，只有这么多。我考虑使用 sample() 进行自助取样，但必须回到我的笔记中，看看为什么我没有采取这种方法...也许这部分讨论应该转到 CrossValidated 上进行。 - avg

问题是你希望从随机数中推导出什么。你的样本可能太小，无法从你的方法中得出任何有意义的结论。 - Roland

Answer 1

你的样本大小有点不确定是否适合混合分布，但是没有关系。你可以按照以下方式从拟合的混合分布中抽样：

probs <- dat.demand2050.mixmdl$lambda
m <- dat.demand2050.mixmdl$mu
s <- at.demand2050.mixmdl$sigma

N <- 1e5
grp <- sample(length(probs), N, replace=TRUE, prob=probs)
x <- rnorm(N, m[grp], s[grp])

- Hong Ooi

---

你的方法似乎过分强调了与Roland解决方案相同的低分布。将你的输出密度与起始密度和@CnrL的解决方案的输出进行比较。这段代码看起来没问题，但结果似乎有误。我不确定为什么。 - David Marx

结果与@CnrL的完全相同。使用N = 1e5运行他们的解决方案。至于起始密度，谁知道15个数据点会发生什么。 - Hong Ooi

@DavidMarx 两种解决方案都没有给出与原始样本相同的密度图。这是一个样本大小的问题。 - Roland

请看我在问题上的评论，了解为什么我正在使用小尺寸。我已经使用较大的N（请检查我之前评论中的链接）运行了@CnrL的解决方案，但仍然给出了较低的峰值。 - avg

@HongOoi，我使用您的解决方案得到的结果与CrnL的解决方案不同。以下是我从您的解决方案（N=1e5）、他的解决方案（N=1e5）和原始小数据集中得到的密度的比较：http://i.imgur.com/cMKnhhf.jpg - David Marx

Answer 2

你的方法是正确的。

对于来自混合分布的每个样本，你只需要选择该样本应来自哪个高斯分布组件，然后从该分布中抽取样本。

你可以使用已找到的混合比例来选择两个分布之间：模拟一个介于0和1之间的随机数，如果随机数小于第一个比例，则从第一个分布中进行采样，否则从第二个分布中进行采样。

最后，使用rnorm函数从相关的高斯分布中进行采样。

dat.demand2050.mixmdl.prop=c(0.2783939,0.7216061)
dat.demand2050.mixmdl.means=c(56.21150,73.08389)
dat.demand2050.mixmdl.dev=c(3.098292,6.413906)

sampleMixture=function(prop,means,dev){
    # Generate a uniformly distributed random number between 0 and 1
    # in order to choose between the two component distributions
    distTest=runif(1)
    if(distTest<prop[1]){
        # Then sample from the first component of the mixture
        sample=rnorm(1,mean=means[1],sd=dev[1])
    }else{
        # Sample from the second component of the mixture
        sample=rnorm(1,mean=means[2],sd=dev[2])
    }
    return(sample)
}

# Generate a single sample
sampleMixture(dat.demand2050.mixmdl.prop,dat.demand2050.mixmdl.means,dat.demand2050.mixmdl.dev)

# Generate 100 samples and plot resulting distribution
samples=replicate(100,sampleMixture(dat.demand2050.mixmdl.prop,dat.demand2050.mixmdl.means,dat.demand2050.mixmdl.dev))
plot(density(samples))

- CnrL

---

不用谢。不，它并不意味着那样。这是由于 "if" 条件所造成的。使用runif()只是为了在分布选择中引入随机性的一种方式。runif()返回值小于0.28的概率恰好为28%（反之大于等于0.28的概率为72%）。通过检查 runif 是否大于或小于第一个比例（在这种情况下为0.28）并相应地选择混合的第一或第二个组件，我们正确地加权了概率。 - CnrL

谢谢，你的解决方案似乎很有效。然而，使用 runif() 作为 distTest 的选择是否意味着该值来自两个分布的概率相等，但数据（和拟合）表明“概率”约为0.3和0.7？ - avg

你应该避免使用循环。从两个正态分布和均匀分布中各创建100个样本，并使用ifelse。 - Roland

@AdvaitGodbole 这些概率在 if 语句中被考虑进去了。该函数从均匀分布中进行采样，以便我们随机选择其中一个混合物，但是该选择将按照这些概率指定的方式发生。 - David Marx

@Roland，您上面的建议肯定会使事情更快，下面的答案非常优雅，但是对于这个问题，我更喜欢上面不那么快速的代码，因为我认为它更清楚地向OP解释了采样的工作原理。 - CnrL