在 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. 如果我的方法是正确的，我如何使用这个结果从这个混合分布中生成一组随机值？

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

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