Вопрос

Допустим, у меня неизвестная плотность a .

Все, что я знаю, это сетка вероятностей ( probs ) квантилей ( quants ).

Как я могу генерировать случайные выборки из неизвестной плотности?

Это то, что у меня есть до сих пор.

Я даю отбраковку выборки, но я не привязан к этому методу. Здесь я подгоняю многочлен (6 градусов) к квантилям. Целью этого является преобразование дискретных квантилей в гладкую непрерывную функцию. Это дает мне эмпирический CDF. Затем я использую выборку отклонения, чтобы получить фактические выборки из CDF. Есть ли в R удобный способ преобразования выборок из CDF в выборки плотности, или я сделал это запутанным способом, когда есть лучшая альтернатива?

 # unknown and probably not normal, but I use rnorm here because it is easy
a <- c(exp(rnorm(200, 5, .8)))
probs <- seq(0.05, 0.95, 0.05)
quants <- quantile(a, probs)
df_quants <- tibble::tibble(cum_probs, quants)
df_quants <- df_quants
fit <- lm(quants ~ poly(cum_probs, 6), df_quants)
df_quants$fit <- predict(fit, df_quants)

p <- df_quants %>%
  ggplot(aes(x = cum_probs, y = quants)) 
  geom_line(aes(y = quants), color = "black", size = 1)  
  geom_line(aes(y = fit), color = "red", size = 1)
  

CDF

 count = 1
accept = c()
X <- runif(50000, 0, 1)
U <- runif(50000, 0, 1)
estimate <- function(x){
  new_x <- predict(fit, data.frame(cum_probs = c(x)))
  return(new_x)
while(count <= 50000 amp; length(accept) < 40000){
  test_u = U[count]
  test_x = estimate(X[count])/(1000*dunif(X[count], 0, 1))
  if(test_u <= test_x){
    accept = rbind(accept, X[count])
    count = count   1
  }
    count = count   1
}
p2 <- as_tibble(accept, name = V1) %>%
  ggplot(aes(x = V1))  
  geom_histogram(bins = 45)
}
  

Образцы CDF

Я не думаю, что выборка отклонения необходима, с помощью подгонки Bspline я смог генерировать разумные выборки с помощью обратного преобразования, но мне также нужна сетка с более высоким разрешением. Хвосты немного отклонены.

Предположение, которое я делаю здесь, состоит в том, что подгонка Bspline к плотной сетке квантилей приближает обратную функцию CDF. Как только эта кривая fut, я могу просто использовать случайные формы U[0,1]

 library(splines2)

a <- c(exp(rnorm(200, 5, .8)))
cum_probs <- seq(0.01, 0.99, 0.01)
quants <- quantile(a, cum_probs)
df_quants <- tibble::tibble(cum_probs, quants)
fit_spline <- lm(quants ~ bSpline(cum_probs, df = 9), df_quants)
df_quants$fit_spline <- predict(fit_spline, df_quants)
estimate <- function(x){
  new_x <- predict(fit_spline, data.frame(cum_probs = c(x)))
  return(new_x)
}
e <- runif(10000, 0, 1)
y <-(estimate(e))
df_density <- tibble(y)
df_densitya <- tibble(a)
py <- df_density %>%
  ggplot(aes(x = y))  
  geom_histogram()
pa <- df_densitya %>%
  ggplot(aes(x = a))  
  geom_histogram(bins = 45)
  

исходная плотность

Примеры обратного преобразования

сводная статистика

оригинальный dist a

 Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
20.36   80.84  145.25  195.72  241.22 1285.24
  

генерируется из квантилей y

 Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
28.09   81.78  149.53  189.07  239.62  667.27