2013-06-13 20 views
9

Załóżmy, że mam zestaw liczb, które, jak podejrzewam, pochodzą z tego samego rozkładu.Generowanie losowej liczby z obiektu gęstości (lub szerzej z zestawu liczb)

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

Chciałbym funkcji, która losowo generuje liczbę z tego samego nieznanego rozkładu. Jedna z metod, o której myślałem, to stworzenie obiektu o nazwie density, a następnie pobranie z niego CDF i przyjęcie odwrotnego CDF losowej zmiennej jednolitej (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 

Moje pytania są następujące:

  1. Czy istnieje lepsza (lub wbudowane R) sposób, aby wygenerować liczbę losową z obiektu gęstości?
  2. Czy są jakieś inne pomysły na generowanie losowej liczby z zestawu liczb (oprócz sample)?
+0

Teoria tego jest znacznie bardziej subtelna. Jak szacuje się gęstość? Które jądro jest używane? Czy wokół tego oszacowania są przedziały ufności? Czy to może być model mieszany? itp. –

Odpowiedz

9

Aby wygenerować dane z oszacowania gęstości, wystarczy losowo wybrać jeden z oryginalnych punktów danych i dodać losowy "błąd" na podstawie jądra z oszacowania gęstości, dla domyślnego "Gaussa" oznacza to tylko wybór element losowy z oryginalnego wektora i dodanie losowego normalnego ze średnią 0 i sd równa szerokości pasma stosowanego:

den <- density(x) 

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

Innym rozwiązaniem jest wyposażenie gęstości przy użyciu funkcji logspline pakiecie logspline (stosuje się inną metodę oszacowanie gęstości), a następnie użyj funkcji rlogspline w tym pakiecie, aby wygenerować nowe dane z szacowanej gęstości.

2

Jeśli wszystko, czego potrzebujesz, to narysować wartości z istniejącej puli liczb, to jest to droga, którą trzeba przejść, sample.
Jeśli chcesz pobrać z założonej dystrybucji podstawowej, użyj density i dopasuj ją do swojej przypuszczalnej dystrybucji, aby uzyskać niezbędne współczynniki (średnie, SD, itd.) I użyj odpowiedniej funkcji rozkładu R.

Co więcej, zajrzę do Rozdziału 7.3 ("metoda odrzucania") w Receptach Numerycznych w C, aby uzyskać sposoby "selektywnego" badania według dowolnego rozkładu. Kod jest na tyle prosty, że można go łatwo przetłumaczyć na R. Mój zakład jest taki, że ktoś już to zrobił i post zamieści lepszą odpowiedź.

Powiązane problemy