Jak narysować losowe próbki z nieparametrycznego rozkładu szacunkowego?

14

Mam próbkę 100 punktów, które są ciągłe i jednowymiarowe. Oszacowałem jego nieparametryczną gęstość za pomocą metod jądra. Jak narysować losowe próbki z tego szacunkowego rozkładu?

lovekesh
źródło

Odpowiedzi:

21

Szacunkowa gęstość jądra jest rozkładem mieszanki; dla każdej obserwacji jest jądro. Jeśli jądro ma gęstość skalowaną, prowadzi to do prostego algorytmu próbkowania z oszacowania gęstości jądra:

repeat nsim times:
  sample (with replacement) a random observation from the data
  sample from the kernel, and add the previously sampled random observation

hxjaN.(μ=xja,σ=h)

# Original distribution is exp(rate = 5)
N = 1000
x <- rexp(N, rate = 5)

hist(x, prob = TRUE)
lines(density(x))

# Store the bandwith of the estimated KDE
bw <- density(x)$bw

# Draw from the sample and then from the kernel
means <- sample(x, N, replace = TRUE)
hist(rnorm(N, mean = means, sd = bw), prob = TRUE)

M.

M = 10
hist(rnorm(N * M, mean = x, sd = bw))

Jeśli z jakiegoś powodu nie możesz czerpać z jądra (np. Twoje jądro nie jest gęstością), możesz spróbować z ważnością próbkowania lub MCMC . Na przykład przy użyciu próbkowania ważności:

# Draw from proposal distribution which is normal(mu, sd = 1)
sam <- rnorm(N, mean(x), 1)

# Weight the sample using ratio of target and proposal densities
w <- sapply(sam, function(input) sum(dnorm(input, mean = x, sd = bw)) / 
                                 dnorm(input, mean(x), 1))

# Resample according to the weights to obtain an un-weighted sample
finalSample <- sample(sam, N, replace = TRUE, prob = w)

hist(finalSample, prob = TRUE)

PS Z podziękowaniami dla Glen_b, który przyczynił się do odpowiedzi.

Matteo Fasiolo
źródło
1
Niestety przeszedłem do ważnego próbkowania, a potem zdałem sobie sprawę, że zwykle próbkowanie jest prostsze. Dodałem twoje wstępne wyjaśnienie do moich odpowiedzi. Wielkie dzięki
Matteo Fasiolo,
@ Matteo Fasiolo - Czy masz jakieś odniesienia do artykułu, który mogę zacytować dla tej metody?
Pallavi