2009-09-14 10 views
6

Tengo un objeto densidad dd creado de esta manera:Generar desvía al azar estocásticos de un objeto de densidad con R

x1 <- rnorm(1000) 
x2 <- rnorm(1000, 3, 2) 
x <- rbind(x1, x2) 
dd <- density(x) 
plot(dd) 

que produce esta distribución muy no gaussiana:

alt text http://www.cerebralmastication.com/wp-content/uploads/2009/09/nongaus.png

lo haría en última instancia, le gusta obtener desviaciones aleatorias de esta distribución, similar a la forma en que rnorm se desvía de una distribución normal.

La forma en que trato de descifrar esto es obtener la CDF de mi kernel y luego hacer que me diga la variante si le paso una probabilidad acumulativa (CDF inversa). De esa manera puedo convertir un vector de variables aleatorias uniformes en dibujos de la densidad.

Parece que lo que estoy tratando de hacer debe ser algo básico que otros hayan hecho antes que yo. ¿Hay una manera simple o una simple función para hacer esto? Odio reinventar la rueda.

FWIW Encontré this R Help article pero no puedo entender lo que están haciendo y la salida final no parece producir lo que estoy buscando. Pero podría ser un paso en el camino que simplemente no entiendo.

He considerado ir con un Johnson distribution from the suppdists package pero Johnson no me dará la buena joroba bimodal que tienen mis datos.

+0

más una cuestión de estadísticas de la programación ... –

+0

sé las estadísticas. Quiero implementar el método de estadísticas en un idioma determinado. Eso es programación. –

Respuesta

2

Esto es solo una mezcla de normales. Entonces, ¿por qué no algo como:

rmnorm <- function(n,mean, sd,prob) { 
    nmix <- length(mean) 
    if (length(sd)!=nmix) stop("lengths should be the same.") 
    y <- sample(1:nmix,n,prob=prob, replace=TRUE) 
    mean.mix <- mean[y] 
    sd.mix <- sd[y] 
    rnorm(n,mean.mix,sd.mix) 
} 
plot(density(rmnorm(10000,mean=c(0,3), sd=c(1,2), prob=c(.5,.5)))) 

Esto debería estar bien si todo lo que necesita son muestras de esta distribución de mezcla.

+0

¡Me gusta la idea! Pero mi ejemplo se simplifica demasiado con fines ilustrativos. En realidad, no conozco los dos modos y podría tener solo un modo y una cola larga (es decir, leptokurtosis). Pero me gusta tu ejemplo. No podría haber programado eso tan brevemente. BTW, creo que faltan CA en: diagrama (densidad (rmnorm (10000, mean = c (0,3), sd = c (1,2), prob = c (.5, .5)))) –

+0

@JD Long: gracias por detectar el error tipográfico. –

+1

Por eso quiere la respuesta de Hadley: remuestrearla. Recuerde que su gráfica de densidad es/solo una estimación/que también depende de su parámetro de suavizado. –

9

enfoque alternativo:

sample(x, n, replace = TRUE) 
+1

bootstrapping ftw! –

+0

Sí, he estado pensando en esto. Si hago una muestra + un sorteo de una normal debería terminar engrosando mis colas de la misma manera que lo haría el núcleo, ¿verdad? Suponiendo que param mi normal de la misma manera que lo haría el método kerneling. –

+2

Sí, agregue rvs normales con media cero y sd = ancho de banda de densidad estimada: muestra (x, n, reemplazar = TRUE) + rnorm (n, 0, sd = 0.4214) Simulando así se discute en el libro de Silverman de 1986 estimación de densidad. –

Cuestiones relacionadas