Génération de variables aléatoires à partir d'un mélange de distributions normales

20

Comment puis-je échantillonner à partir d'une distribution de mélange, et en particulier d'un mélange de distributions normales dans R? Par exemple, si je voulais échantillonner à partir de:

0,3×N(0,1)+0,5×N(dix,1)+0,2×N(3,.1)

comment pourrais-je faire ça?

gung - Réintégrer Monica
la source
3
Je n'aime vraiment pas cette façon de désigner un mélange. Je sais que c'est conventionnellement fait comme ça, mais je trouve cela trompeur.La notation suggère que pour échantillonner, vous devez échantillonner les trois normales et peser les résultats par ces coefficients qui ne seraient évidemment pas corrects. Quelqu'un connaît une meilleure notation?
StijnDeVuyst
Je n'ai jamais eu cette impression. Je pense aux distributions (dans ce cas les trois distributions normales) comme des fonctions et le résultat est une autre fonction.
roundsquare
@StijnDeVuyst vous voudrez peut-être visiter cette question provient de votre commentaire: stats.stackexchange.com/questions/431171/…
ankii
@ankii: merci de l'avoir signalé!
StijnDeVuyst

Réponses:

32

Il est recommandé d'éviter les forboucles Rpour des raisons de performances. Une solution alternative qui exploite le fait rnormest vectorisée:

N <- 100000

components <- sample(1:3,prob=c(0.3,0.5,0.2),size=N,replace=TRUE)
mus <- c(0,10,3)
sds <- sqrt(c(1,1,0.1))

samples <- rnorm(n=N,mean=mus[components],sd=sds[components])
M. Berk
la source
3
Vous pouvez également utiliser les propriétés de la distribution normale pour remplacer la dernière ligne par samples <- rnorm(N)*sds[components]+mus[components]. Je trouve plus facile à lire :)
Elvis
Très élégant (cc @Elvis)!
Itamar
18

En général, l'un des moyens les plus simples d'échantillonner à partir d'une distribution de mélange est le suivant:

Étapes de l'algorithme

1) Générer une variable aléatoire UUniforme(0,1)

U[je=1kpk,je=1k+1pk+1)pkkthkth

3) Répétez les étapes 1) et 2) jusqu'à ce que vous ayez la quantité désirée d'échantillons de la distribution du mélange

Maintenant, en utilisant l'algorithme général donné ci-dessus, vous pouvez échantillonner à partir de votre exemple de mélange de normales en utilisant le Rcode suivant :

#The number of samples from the mixture distribution
N = 100000                 

#Sample N random uniforms U
U =runif(N)

#Variable to store the samples from the mixture distribution                                             
rand.samples = rep(NA,N)

#Sampling from the mixture
for(i in 1:N){
    if(U[i]<.3){
        rand.samples[i] = rnorm(1,0,1)
    }else if(U[i]<.8){
        rand.samples[i] = rnorm(1,10,1)
    }else{
        rand.samples[i] = rnorm(1,3,.1)
    }
}

#Density plot of the random samples
plot(density(rand.samples),main="Density Estimate of the Mixture Model")

#Plotting the true density as a sanity check
x = seq(-20,20,.1)
truth = .3*dnorm(x,0,1) + .5*dnorm(x,10,1) + .2*dnorm(x,3,.1)
plot(density(rand.samples),main="Density Estimate of the Mixture Model",ylim=c(0,.2),lwd=2)
lines(x,truth,col="red",lwd=2)

legend("topleft",c("True Density","Estimated Density"),col=c("red","black"),lwd=2)

Ce qui génère:

entrez la description de l'image ici

et comme vérification de santé mentale:

entrez la description de l'image ici


la source
Salut! Merci beaucoup! Cette réponse m'a beaucoup aidé. Je l'utilise dans un projet de recherche. Je souhaite citer une référence pour ce qui précède. Pouvez-vous s'il vous plaît suggérer une citation d'article de recherche.
Abhishek Bhatia
7

kR

set.seed(8)               # this makes the example reproducible
N     = 1000              # this is how many data you want
probs = c(.3,.8)          # these are *cumulative* probabilities; since they 
                          #   necessarily sum to 1, the last would be redundant
dists = runif(N)          # here I'm generating random variates from a uniform
                          #   to select the relevant distribution

# this is where the actual data are generated, it's just some if->then
#   statements, followed by the normal distributions you were interested in
data = vector(length=N)
for(i in 1:N){
  if(dists[i]<probs[1]){
    data[i] = rnorm(1, mean=0, sd=1)
  } else if(dists[i]<probs[2]){
    data[i] = rnorm(1, mean=10, sd=1)
  } else {
    data[i] = rnorm(1, mean=3, sd=.1)
  }
}

# here are a couple of ways of looking at the results
summary(data)
#    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# -3.2820  0.8443  3.1910  5.5350 10.0700 13.1600 

plot(density(data))

entrez la description de l'image ici

gung - Réintégrer Monica
la source
Belle réponse, vous m'avez battu pour poster: P
1
Merci pour l'astuce, @BabakP. Je ne sais pas ce que c'était. C'était quelque chose dans la ifelse()déclaration, mais je devrai le découvrir plus tard. J'ai remplacé ce code avec une boucle.
gung - Rétablir Monica
6
RfindInterval()cumsum()μmuσ2spmix <- function(n,mu,s,p) { ii <- findInterval(runif(n),cumsum(p))+1; x <- rnorm(n,mean=mu[ii],sd=sqrt(s[ii])); return(x); }
1
@Macro, code très vrai et très agréable! Je n'ai jamais vu la findInterval()commande auparavant, cependant, j'aime écrire du code ici aussi simplement que possible car je veux que ce soit un outil de compréhension plutôt que d'efficacité.
1
J'ai dit que c'étaient de bonnes réponses. Mon but n'était pas de vous critiquer mais de proposer une approche qui se généralise facilement à plus de trois dimensions en ne modifiant qu'un seul argument, pas n'importe quel code. Je ne comprends pas pourquoi ce que vous avez écrit est plus transparent que ce que j'ai écrit, mais je ne veux certainement pas contester cela. À votre santé.
Macro
0

Déjà donné des réponses parfaites, donc pour ceux qui veulent y parvenir en Python, voici ma solution:

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

mu = [0, 10, 3]
sigma = [1, 1, 1]
p_i = [0.3, 0.5, 0.2]
n = 10000

x = []
for i in range(n):
    z_i = np.argmax(np.random.multinomial(1, p_i))
    x_i = np.random.normal(mu[z_i], sigma[z_i])
    x.append(x_i)

def univariate_normal(x, mean, variance):
    """pdf of the univariate normal distribution."""
    return ((1. / np.sqrt(2 * np.pi * variance)) * 
            np.exp(-(x - mean)**2 / (2 * variance)))

a = np.arange(-7, 18, 0.01)
y = p_i[0] * univariate_normal(a, mean=mu[0], variance=sigma[0]**2) + p_i[1] * univariate_normal(a, mean=mu[1], variance=sigma[0]**2)+ p_i[2] * univariate_normal(a, mean=mu[2], variance=sigma[0]**2)

fig, ax = plt.subplots(figsize=(8, 4))

ax.hist(x, bins=100, density=True)
ax.plot(a, y)

entrez la description de l'image ici

UN RAT
la source