Generare un numero casuale da un oggetto densità (o più genericamente da un insieme di numeri)

Question

Diciamo che ho una serie di numeri che sospetto provengano dalla stessa distribuzione.

set.seed(20130613) 
x <- rcauchy(10) 

Vorrei una funzione che generi a caso un numero dalla stessa distribuzione sconosciuta. Un approccio che ho pensato è quello di creare un oggetto density e quindi ottenere il CDF da quello e prendere il CDF inverso di una variabile uniforme casuale (see Wikipedia).

den <- density(x) 

#' Generate n random numbers from density() object 
#' 
#' @param n The total random numbers to generate 
#' @param den The density object from which to generate random numbers 
rden <- function(n, den) 
{ 
     diffs <- diff(den$x) 
     # Making sure we have equal increments 
     stopifnot(all(abs(diff(den$x) - mean(diff(den$x))) < 1e-9)) 
     total <- sum(den$y) 
     den$y <- den$y/total 
     ydistr <- cumsum(den$y) 
     yunif <- runif(n) 
     indices <- sapply(yunif, function(y) min(which(ydistr > y))) 
     x <- den$x[indices] 

     return(x) 
} 

rden(1, den) 
## [1] -0.1854121 

Le mie domande sono le seguenti:

1. C'è un modo migliore (o integrato in R) per generare un numero casuale da un oggetto densità?
2. Ci sono altre idee su come generare un numero casuale da un insieme di numeri (oltre a sample)?

Answer 1

Per generare dati da una stima della densità, basta scegliere casualmente uno dei punti dati originali e aggiungere un pezzo di "errore" casuale basato sul kernel dalla stima della densità, per il valore predefinito di "Gaussian" questo significa semplicemente scegliere un elemento casuale dal vettore originale e aggiungere un casuale normale con media 0 e sD uguale alla larghezza di banda utilizzata:

den <- density(x) 

N <- 1000 
newx <- sample(x, N, replace=TRUE) + rnorm(N, 0, den$bw) 

Un'altra opzione è quella di montare una densità utilizzando la funzione logspline dal pacchetto logspline (utilizza un diverso metodo di stimare una densità), quindi utilizzare la funzione rlogspline in quel pacchetto per generare nuovi dati dalla densità stimata.

Answer 2

Se tutto ciò che serve è disegnare valori dal pool esistente di numeri, quindi sample è la strada da percorrere.
Se si desidera disegnare dalla presunta distribuzione sottostante, utilizzare density e adattarlo alla distribuzione presunta per ottenere i coefficienti necessari (mean, sd, ecc.) E utilizzare la funzione di distribuzione R appropriata.

Oltre a ciò, darei un'occhiata al Capitolo7.3 ("metodo di rifiuto") di Ricette numeriche in C per modi di campionare "selettivamente" in base a qualsiasi distribuzione. Il codice è abbastanza semplice da essere facilmente tradotto in R. La mia scommessa è già stata fatta da qualcuno e pubblicheremo una risposta migliore di questa.