R : функция для создания распределения смеси

Мне нужно создать образцы из смешанного дистрибутива

  • 40% выборок исходят из гауссова (среднее значение = 2, стандартное отклонение = 8)

  • 20% образцов поступают из Коши (местоположение = 25, масштаб = 2)

  • 40% выборок исходят от гауссова (среднее = 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)

Я получаю подобную гистограмму (которая, как я знаю, неверна) - гистограмма предполагаемого распределения

Что я делаю не так? Может ли кто-нибудь дать несколько указателей, пожалуйста?


r probability mcmc markov-chains mixture-model
person Raaj    schedule 05.05.2014    source источник
comment
только одно я заметил: вы можете создать тот же сюжет, используя hist(dmix(seq(-5,5,by = 0.01)))без зацикливания и vector   -  person talat    schedule 06.05.2014
comment
Вы можете сделать гистограмму случайных выборок? 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()   -  person rawr    schedule 06.05.2014
comment
@rawr Я мог бы использовать это - dmix2 ‹- function(x){ prob ‹- c(rnorm(x * .4, среднее = 2, sd = 8), rcauchy(x * .2, location = 25, scale = 2) ),rnorm(x * .4, mean = 10, sd = 6)) return (prob) } hist(prob) Но мне нужно получить вероятность того, что выборка принадлежит распределению. dnorm, dcauchy и т. д. дают вероятность. Я полагаю, что тогда мне следует использовать rnorm и rcauchy. [Здесь нельзя добавить график гистограммы]   -  person Raaj    schedule 06.05.2014

Ответы (2)


arrow_upward
12
arrow_downward

Конечно, есть и другие способы сделать это, но пакет distr делает это чертовски простым. (См. также этот ответ для другого примера и более подробной информации о 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")

введите здесь описание изображения

person Josh O'Brien    schedule 05.05.2014
comment
Кажется, это то, что я ищу. - person Raaj; 06.05.2014
comment
О, это именно то, что я ищу! Большое спасибо! Не знал о существовании дистрибутива. Наткнулся на микстулы. Однако это, кажется, больше для анализа, чем для генерации. - person Raaj; 06.05.2014
comment
@Радж Отлично. Рад, что помог. Как я уверен, вы только что видели, это лишь малая часть того, что Ruckdeschel et al. собрались! - person Josh O'Brien; 06.05.2014
comment
Могу я вас побеспокоить еще одним вопросом? Мне кажется, что [-100 to 100] похоже на усечение (или это одно и то же?) Какова его цель? Я, конечно, изучу документацию, но, если можно, дайте мне знать. Спасибо. - person Raaj; 06.05.2014
comment
когда вы рисуете случайные величины Коши, если вы не обрезаете хвосты таким образом, вы получаете гистограмму со смехотворно большим x-диапазоном (например, от -1e-7 до +1e7), и тогда вы не можете видеть тело дистрибутива очень хорошо вообще. - person Ben Bolker; 06.05.2014
comment
А, спасибо. Это ускользнуло от меня на мгновение. Поиграем с параметрами и понаблюдаем еще раз, чтобы отпечататься. - person Raaj; 06.05.2014

arrow_upward
2
arrow_downward

Всегда используйте 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
person Meme    schedule 27.10.2016