2009-09-14 9 views
6

J'ai un objet de densité dd créé comme ceci:Générer déviants aléatoires stochastiques à partir d'un objet de densité avec R

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

qui produit cette distribution très non gaussienne:

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

Je voudrais En fin de compte, il est préférable d'obtenir des déviations aléatoires de cette distribution, de la même façon que rnorm s'écarte d'une distribution normale.

La façon dont j'essaie de résoudre ceci est d'obtenir le CDF de mon noyau et de l'obtenir pour me dire la variance si je lui transmets une probabilité cumulée (CDF inverse). De cette façon, je peux transformer un vecteur de variables aléatoires uniformes en tirages à partir de la densité.

Il semble que ce que j'essaie de faire devrait être quelque chose de fondamental que les autres ont fait avant moi. Y a-t-il un moyen simple ou une fonction simple de le faire? Je déteste réinventer la roue. FWIW J'ai trouvé this R Help article mais je ne peux pas dire ce qu'ils font et la sortie finale ne semble pas produire ce que je cherche. Mais cela pourrait être une étape que je ne comprends tout simplement pas.

J'ai envisagé d'aller avec un Johnson distribution from the suppdists package mais Johnson ne me donnera pas la bosse bimodale sympa que mes données ont.

+0

plus d'une question que la programmation des statistiques ... –

+0

Je sais que les stats. Je veux implémenter la méthode de statistiques dans une langue donnée. C'est la programmation. –

Répondre

2

Ceci est juste un mélange de normales. Alors pourquoi pas quelque chose comme:

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

Cela devrait aller si tout ce dont vous avez besoin sont des échantillons de cette distribution de mélange.

+0

J'aime l'idée! Mais mon exemple est trop simplifié à des fins d'illustration. En réalité, je ne connais pas les deux modes et il pourrait avoir un seul mode et une longue queue de cul (c'est-à-dire leptokurtosis). Mais j'aime ton exemple. Je n'aurais pas pu programmer cela aussi succinctement. BTW, je crois qu'il manque ac en: plot (densité (rmnorm (10000, moyenne = c (0,3), sd = c (1,2), prob = c (0,5, 0,5)))) –

+0

@JD Long: merci de repérer la faute de frappe. –

+1

C'est pourquoi vous voulez la réponse de Hadley - le rééchantillonnage est. Rappelez-vous que votre diagramme de densité est/juste une estimation/qui dépend également de votre paramètre de lissage. –

9

approche alternative:

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

bootstrapping ftw! –

+0

ouais, j'ai trop pensé ça. Si je fais un échantillon + un tirage d'une normale je devrais finir par épaissir mes queues de la même façon que le noyau, non? En supposant que je param ma normale de la même manière la méthode de kerneling serait. –

+2

Oui, ajoutez des valeurs normales avec zéro moyen et sd = bande passante à partir de l'estimation de densité: échantillon (x, n, replace = TRUE) + rnorm (n, 0, sd = 0,4214) La simulation comme ceci est discutée dans le livre de Silverman de 1986 sur estimation de densité. –