Jak najlepiej symulować dowolną losową zmienną jednowymiarową za pomocą jej funkcji prawdopodobieństwa?

Question

W języku R, jaki jest najlepszy sposób symulacji arbitralnej zmiennej losowej jednowymiarowej, jeśli dostępna jest tylko funkcja gęstości prawdopodobieństwa?

Answer 1

Oto (powolna) implementacja metody odwrotnego cdf, gdy otrzymuje się tylko gęstość.

den<-dnorm #replace with your own density 

#calculates the cdf by numerical integration 
cdf<-function(x) integrate(den,-Inf,x)[[1]] 

#inverts the cdf 
inverse.cdf<-function(x,cdf,starting.value=0){ 
lower.found<-FALSE 
lower<-starting.value 
while(!lower.found){ 
    if(cdf(lower)>=(x-.000001)) 
    lower<-lower-(lower-starting.value)^2-1 
    else 
    lower.found<-TRUE 
} 
upper.found<-FALSE 
upper<-starting.value 
while(!upper.found){ 
    if(cdf(upper)<=(x+.000001)) 
    upper<-upper+(upper-starting.value)^2+1 
    else 
    upper.found<-TRUE 
} 
uniroot(function(y) cdf(y)-x,c(lower,upper))$root 
} 

#generates 1000 random variables of distribution 'den' 
vars<-apply(matrix(runif(1000)),1,function(x) inverse.cdf(x,cdf)) 
hist(vars) 

Answer 2

Zastosowanie skumulowany rozkład funkcji http://en.wikipedia.org/wiki/Cumulative_distribution_function

Następnie wystarczy użyć jego odwrotność. Sprawdź tutaj lepszy obraz http://en.wikipedia.org/wiki/Normal_distribution

To znaczy: podnieść liczbę losową z [0,1] i ustawić jako CDF, a następnie sprawdzić Wartość

nazywany jest również odwrotna dystrybuanta.

Answer 3

Możesz użyć metropolis-hastings, aby pobrać próbki z gęstości.

Answer 4

Aby wyjaśnić "użytkowania Metropolis-Hastings" odpowiedź powyżej:

Załóżmy ddist() jest twoja funkcja gęstości prawdopodobieństwa

coś takiego:

n <- 10000 
cand.sd <- 0.1 
init <- 0 
vals <- numeric(n) 
vals[1] <- init 
oldprob <- 0 
for (i in 2:n) { 
    newval <- rnorm(1,mean=vals[i-1],sd=cand.sd) 
    newprob <- ddist(newval) 
    if (runif(1)<newprob/oldprob) { 
     vals[i] <- newval 
    } else vals[i] <- vals[i-1] 
    oldprob <- newprob 
} 

Uwagi:

1. całkowicie nieprzetestowana
2. wydajność zależy od dystrybucji kandydata (tj. wartość cand.sd). Dla maksymalnej wydajności, dostroić cand.sd do stopy akceptacji 25-40%
3. wyniki zostaną autokorelacji ... (chociaż myślę, że zawsze można sample() wyniki ich wspiąć lub cienkiej)
4. mogą wymagać odrzucić „wypalenia”, jeśli wartość początkowa jest dziwne

klasycznym podejściem do tego problemu jest odrzucenie próbkowania (patrz np Press i wsp Numerical Recipes)