Menghasilkan variabel acak dari campuran distribusi Normal

20

Bagaimana saya bisa mengambil sampel dari distribusi campuran, dan khususnya campuran distribusi normal R? Misalnya, jika saya ingin mengambil sampel dari:

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

bagaimana saya bisa melakukan itu?

gung - Reinstate Monica
sumber
3
Saya benar-benar tidak suka cara ini menunjukkan campuran. Saya tahu ini dilakukan secara konvensional seperti ini, tetapi saya merasa itu menyesatkan. Notasi menunjukkan bahwa untuk sampel, Anda perlu sampel ketiga normals dan menimbang hasilnya dengan koefisien-koefisien yang jelas tidak benar. Adakah yang tahu notasi yang lebih baik?
StijnDeVuyst
Saya tidak pernah mendapat kesan seperti itu. Saya memikirkan distribusi (dalam hal ini tiga distribusi normal) sebagai fungsi dan kemudian hasilnya adalah fungsi lain.
roundsquare
@StijnDeVuyst Anda mungkin ingin mengunjungi pertanyaan ini berasal dari komentar Anda: stats.stackexchange.com/questions/431171/…
ankii
@ankii: terima kasih sudah menunjukkannya!
StijnDeVuyst

Jawaban:

32

Ini praktik yang baik untuk menghindari forloop Rkarena alasan kinerja. Solusi alternatif yang mengeksploitasi fakta rnormadalah vektor:

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
sumber
3
Atau, Anda dapat menggunakan properti dari distribusi normal untuk mengganti baris terakhir dengan samples <- rnorm(N)*sds[components]+mus[components]. Saya merasa lebih mudah untuk membaca :)
Elvis
Sangat elegan (cc @ Elvis)!
Itamar
18

Secara umum, salah satu cara termudah untuk mengambil sampel dari distribusi campuran adalah sebagai berikut:

Langkah-langkah Algoritma

1) Hasilkan variabel acak USeragam(0,1)

U[saya=1khalk,saya=1k+1halk+1)halkkthkth

3) Ulangi langkah 1) dan 2) hingga Anda memiliki jumlah sampel yang diinginkan dari distribusi campuran

Sekarang menggunakan algoritma umum yang diberikan di atas, Anda dapat mengambil sampel dari contoh campuran normals Anda dengan menggunakan Rkode berikut :

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

Yang menghasilkan:

masukkan deskripsi gambar di sini

dan sebagai cek kewarasan:

masukkan deskripsi gambar di sini


sumber
Hai! Terima kasih banyak! Jawaban ini sangat membantu saya. Saya menggunakan ini dalam proyek penelitian. Saya ingin mengutip referensi untuk hal di atas. Bisakah Anda menyarankan kutipan artikel penelitian.
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))

masukkan deskripsi gambar di sini

gung - Reinstate Monica
sumber
Jawaban yang bagus, Anda mengalahkan saya untuk memposting: P
1
Terima kasih atas tipnya, @BabakP. Saya tidak yakin apa itu. Itu adalah sesuatu dalam ifelse()pernyataan itu, tetapi saya harus mencari tahu nanti. Saya mengganti kode itu dengan loop.
gung - Reinstate 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, kode yang sangat benar dan sangat bagus! Saya belum pernah melihat findInterval()perintah sebelumnya, namun, saya suka menulis kode di sini sesederhana mungkin karena saya ingin itu menjadi alat untuk memahami daripada efisiensi.
1
Saya mengatakan ini adalah jawaban yang bagus. Tujuan saya bukan untuk mengkritik Anda tetapi untuk menawarkan pendekatan yang dengan mudah digeneralisasikan ke lebih dari tiga dimensi dengan hanya mengubah satu argumen, bukan kode apa pun. Tidak jelas bagi saya mengapa apa yang Anda tulis lebih transparan daripada apa yang saya tulis tetapi saya tentu tidak ingin berdebat tentang itu. Bersulang.
Makro
0

Sudah memberikan jawaban sempurna, jadi bagi mereka yang ingin mencapai ini dengan Python, inilah solusi saya:

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)

masukkan deskripsi gambar di sini

ARAT
sumber