Установка мультимодальных распределений в 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))