Generieren von Zufallsvariablen aus einer Mischung von Normalverteilungen

20

Wie kann ich aus einer Mischungsverteilung und insbesondere einer Mischung von Normalverteilungen in probieren R? Zum Beispiel, wenn ich probieren wollte aus:

0.3×N(0,1)+0.5×N(10,1)+0.2×N(3,.1)

wie könnte ich das machen

gung - Wiedereinsetzung von Monica
quelle
3
Ich mag diese Art, eine Mischung zu bezeichnen, wirklich nicht. Ich weiß, dass es konventionell so gemacht wird, aber ich finde es irreführend. Die Notation legt nahe, dass Sie zum Abtasten alle drei Normalen abtasten und die Ergebnisse mit denjenigen Koeffizienten abwägen müssen, die offensichtlich nicht korrekt wären. Kennt jemand eine bessere Notation?
StijnDeVuyst
Ich habe diesen Eindruck nie bekommen. Ich betrachte die Verteilungen (in diesem Fall die drei Normalverteilungen) als Funktionen und dann ist das Ergebnis eine andere Funktion.
roundsquare
@StijnDeVuyst Vielleicht möchten Sie diese Frage besuchen entstand aus Ihrem Kommentar: stats.stackexchange.com/questions/431171/…
ankii
@ankii: danke für den Hinweis!
StijnDeVuyst

Antworten:

32

Es ist forempfehlenswert, Raus Leistungsgründen keine Schleifen zu bilden. Eine alternative Lösung, die diese Tatsache ausnutzt, rnormist vektorisiert:

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
quelle
3
Alternativ können Sie die Eigenschaften der Normalverteilung verwenden, um die letzte Zeile durch zu ersetzen samples <- rnorm(N)*sds[components]+mus[components]. Ich finde es einfacher zu lesen :)
Elvis
Sehr elegant (cc @ Elvis)!
Itamar
18

Im Allgemeinen ist einer der einfachsten Wege, eine Probe aus einer Gemischverteilung zu entnehmen, der folgende:

Algorithmusschritte

1) Erzeuge eine ZufallsvariableUUniform(0,1)

2) Wenn Intervall, in dem der Wahrscheinlichkeit der -Komponente des Mischungsmodells entsprechen, dann aus der Verteilung der -Komponente erzeugenU[ich=1kpk,ich=1k+1pk+1)pkkthkth

3) Wiederholen Sie die Schritte 1) und 2), bis Sie die gewünschte Menge an Proben aus der Gemischverteilung erhalten haben

Unter Verwendung des oben angegebenen allgemeinen Algorithmus können Sie aus Ihrer Beispielmischung von Normalen mit dem folgenden RCode ein Beispiel erstellen:

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

Welches erzeugt:

Bildbeschreibung hier eingeben

und zur Überprüfung der geistigen Gesundheit:

Bildbeschreibung hier eingeben


quelle
Hallo! Vielen Dank! Diese Antwort hat mir sehr geholfen. Ich benutze dies in einem Forschungsprojekt. Ich möchte eine Referenz für die oben zitieren. Können Sie mir bitte einen Forschungsartikel vorschlagen?
Abhishek Bhatia
7

Konzeptionell wählen Sie mit einiger Wahrscheinlichkeit nur eine Verteilung (aus Möglichkeiten) aus und generieren dann Pseudozufallsvariablen aus dieser Verteilung. In wäre dies (zB): 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))

Bildbeschreibung hier eingeben

gung - Wiedereinsetzung von Monica
quelle
Schöne Antwort, du hast mich zum Posten geschlagen: P
1
Danke für den Tipp, @BabakP. Ich bin mir nicht sicher, was es war. Es war etwas in der ifelse()Aussage, aber ich muss es später herausfinden. Ich habe diesen Code mit einer Schleife ersetzt.
gung - Wiedereinsetzung von Monica
6
(cc @BabakP) Dies sind beide gute Antworten und offensichtlich richtig (+1). Nur ein RProgrammiertrick: Sie können auch die Befehle findInterval()und verwenden cumsum(), um den Code zu vereinfachen und, was noch wichtiger ist, die Verallgemeinerung auf eine andere Anzahl von Dimensionen zu vereinfachen. Zum Beispiel wäre für einen Eingabevektor aus Mittelwerten ( ) und Varianzen ( ) und Mischungswahrscheinlichkeiten ( ) eine einfache Funktion zum Erzeugen von n Abtastwerten aus dieser Mischungμ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); }
Makro
1
@Macro, sehr wahrer und sehr netter Code! Ich habe den findInterval()Befehl vorher noch nicht gesehen , aber ich schreibe hier gerne Code so einfach wie möglich, weil ich möchte, dass er eher zum Verständnis als zur Effizienz beiträgt.
1
Ich sagte, das waren gute Antworten. Mein Ziel war es nicht, Sie zu kritisieren, sondern einen Ansatz anzubieten, der sich leicht auf mehr als drei Dimensionen verallgemeinern lässt, indem nur ein einziges Argument und kein Code geändert wird. Mir ist nicht klar, warum das, was Sie geschrieben haben, transparenter ist als das, was ich geschrieben habe, aber darüber möchte ich sicherlich nicht streiten. Prost.
Makro
0

Bereits perfekte Antworten gegeben, also für diejenigen, die dies in Python erreichen wollen, ist hier meine Lösung:

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)

Bildbeschreibung hier eingeben

EINE RATTE
quelle