Wybrane metody generowania zmiennych losowych

Poniżej przedstawimy algorytmy losowania zmiennych losowych z różnych popularnych i przydatnych rozkładów, których nie omówiliśmy w poprzednim rozdziale.

Proces Poissona

Jednorodny proces Poissona {N(t),t0}\{N(t), t \geq 0\} z częstością λ\lambda, to proces zliczający z niezależnymi przyrostami, spełniający warunki: N(0)=0,P(N(s+t)N(s)=n)=eλt(λt)nn!,     n0,t,s>0. \begin{array}{rcl} N(0)&=&0,\\ P\left(N(s+t) -N(s) = n\right)& = & \frac{e^{\lambda t}(\lambda t)^n}{n!},\ \ \ \ \ n \geq 0, t, s >0.\\ \end{array} Czas oczekiwania TiT_i pomiędzy kolejnymi skokami procesu Poissona ma rozkład wykładniczy E(λ) E(\lambda).

Do generowania trajektorii jednorodnego procesu Poissona często wykorzystuje się jeden z poniższych sposobów. Każdy wykorzystuje inną właściwość procesu Poissona.

Sposób 1

  • Losujemy wektor czasów oczekiwania na kolejne zdarzenie TiT_i z rozkładu wykładniczego.
  • Wyznaczamy czasy pojawiania się kolejnych skoków procesu, jako kumulowane sumy TiT_i Sk=i=1kTi.S_k = \sum_{i=1}^k T_i.
  • Wyznaczamy liczbę skoków do czasu tt N(t)=#{k:Skt}. N(t) = \#\{k : S_k \leq t\}.

Sposób 2

  • Losujemy liczbę skoków procesu do czasu tt, liczba skoków ma rozkład Poissona z parametrem λt\lambda t N(t)=Poiss(λt). N(t) = Poiss(\lambda t).
  • Na odcinku [0,t][0,t] losujemy N(t)N(t) skoków z rozkładu jednostajnego Si=U[0,t],        iN(t). S'_i = U[0,t], \ \ \ \ \ \ \ \ i\leq N(t).

Drugi sposób jest szybszy, ponieważ po pierwszym losowaniu wiemy już ile wartości z rozkładu jednostajnego musimy wylosować. Losowanie z rozkładu jednostajnego jest też szybsze niż losowanie z rozkładu wykładniczego. Poniżej przedstawimy przykładową implementację drugiego sposobu.

rPoisProcess <- function (T=1, lambda=1) {
    N <- rpois(1, T*lambda)
    sort(runif(N,0,T))
}

Przykładowe wywołanie powyższej funkcji.

set.seed(1313)

rPoisProcess(T=10, lambda=0.2)
## [1] 0.6637362

W podobny sposób możemy generować zmienne z jednorodnego pola losowego, czyli z dwuwymiarowego procesu Poissona.

Polem losowym Poissona nazwiemy proces zliczający NN z niezależnymi przyrostami, określony na płaszczyźnie, spełniający warunek: P(N(s1+t1,s2+t2)N(s1,s2)=n)=eλt1t2(λt1t2)nn!. P\left(N(s_1+t_1,s_2+t_2) -N(s_1,s_2) = n\right) = \frac{e^{\lambda t_1 t_2}(\lambda t_1 t_2)^n}{n!}. Liczba zdarzeń w prostokącie o polu SS ma rozkład Poissona o częstości λS\lambda S. Punkty skoków takiego procesu możemy wyznaczyć następująco

  • Losujemy liczbę skoków procesu na prostokącie [0,t1]×[0,t2][0,t_1] \times [0, t_2] N(t1,t2)=Poiss(λt1t2). N(t_1, t_2) = Poiss(\lambda t_1 t_2).
  • Na prostokącie [0,t1]×[0,t2][0,t_1] \times [0, t_2] generujemy N(t1,t2)N(t_1, t_2) punktów z rozkładu jednostajnego.

I implementacja tego algorytmu w programie R.

rRandomField <- function (T1=1, T2=1, lambda=1) {
     N = rpois(1, T1*T2*lambda)
     data.frame(x=runif(N,0,T1), y=runif(N,0,T2))
}

Przykładowe wywołanie powyższej funkcji.

rRandomField(T1=10, T2=10, lambda=0.02)
##          x        y
## 1 5.297950 9.469024
## 2 9.142406 5.634131
## 3 5.115911 6.146278

Rozważmy teraz proces Poissona, którego częstość skoków nie jest stała, ale wynosi λ(t)\lambda(t). Taki proces nazywamy niejednorodnym procesem Poissona. Jest to bardzo często wykorzystywany proces w modelowaniu liczby szkód, częstość występowania szkód najczęściej zmienia się z czasem (liczba oszustw bankomatowych zależy od liczby bankomatów i najczęściej rośnie z czasem, liczba szkód OC zależy od liczby ubezpieczonych samochodów i od szkodowości, oba parametry zmieniają się w czasie itp.).

Załóżmy, że istnieje jakieś maksymalne λmax\lambda_{max}, takie że tλ(t)λmax\forall_t \lambda(t) \leq \lambda_{max}. W takim przypadku trajektorie procesu jednorodnego można przetransformować w trajektorie procesu niejednorodnego.

Taką transformację można wykonać na dwa sposoby.

Sposób 1

  • Wygenerować jednorodny proces Poissona N1N_1 o częstości λmax\lambda_{max}.
  • Każdy punkt SiS_i skoku procesu N1N_1 z prawdopodobieństwem λ(Si)/λmax\lambda(S_i)/\lambda_{max} uznać za punkt skoku procesu N2N_2. Czyli dla każdego punktu z wyjściowego procesu losuje się z prawdopodobieństwem λ(Si)/λmax\lambda(S_i)/\lambda_{max} czy pozostawić ten punkt w nowym procesie. Metoda ta nosi nazwę metody odsiewania. Odrzucając/odsiewając niektóre punkty skoków redukujemy częstość procesu jednorodnego do niejednorodnego.

Sposób 2

  • Wygenerować jednorodny proces Poissona N1N_1 o częstości λ0\lambda_0.
  • Przekształcić czasy skoków SkS_k procesu N1N_1 na czasy SkS_k' procesu N2N_2 tak by 0Skλ(t)dt=0Skλ0dt,\int_0^{S'_k} \lambda(t)dt = \int_0^{S_k} \lambda_0 dt, czyli Sk=Λ1(λ0Sk)S'_k = \Lambda^{-1}(\lambda_0 S_k), gdzie Λ(x)=0xλ(t)dt\Lambda(x) = \int_0^x \lambda(t)dt. Ta metoda nazywana jest metodą zmiany czasu.

Wielowymiarowy rozkład normalny

Gęstość wielowymiarowego rozkładu normalnego Nd(μ,Σ) N_d(\mu, \Sigma) jest postaci f(x)=1(2π)k/2Σ1/2exp(12(xμ)Σ1(xμ)), f(x) = \frac{1}{ (2\pi)^{k/2}|\Sigma|^{1/2} } \exp\Big( {-\tfrac{1}{2}}(x-\mu)'\Sigma^{-1}(x-\mu) \Big), gdzie Σ\Sigma to macierz kowariancji, symetryczna dodatnio określona o wymiarach d×dd \times d a μ\mu to wektor średnich o długości dd.

Zmienne z rozkładu wielowymiarowego normalnego generuje się w dwóch krokach

  • W pierwszym kroku generuje się dd niezależnych zmiennych o rozkładzie normalnym,
  • W drugim kroku przekształca się te zmienne by otrzymać średnią μ\mu i macierz kowariancji Σ\Sigma. Korzysta się z faktu, że jeżeli ZN(0,I)Z\sim N(0, I), to CZ+μ=Nd(μ,CCT). C Z + \mu = N_d(\mu, CC^T).

Wybrane charakterystyki rozkładu ciągłej zmiennej losowej, na przykładzie rozkładu B(2,5)

W takim razie aby wylosować pożądane zmienne wystarczy wylosować wektor niezależnych zmiennych ZZ z następnie przekształcić je liniowo X=CZ+μX = C Z + \mu. W tym celu macierz Σ\Sigma musimy rozłożyć na iloczyn CCTCC^T. Można to zrobić używając różnych faktoryzacji macierzy Σ\Sigma. Wymieńmy trzy

  • Dekompozycja spektralna C=Σ1/2=PΓ1/2P1, C = \Sigma^{1/2} = P \Gamma^{1/2}P^{-1}, gdzie Γ\Gamma to macierz diagonalna wartości własnych na przekątnej a PP to macierz wektorów własnych P1=PTP^{-1}=P^T. Macierze PP i Γ\Gamma dla macierzy Σ\Sigma można znaleźć z użyciem funkcji eigen().

  • Dekompozycja SVD, czyli uogólnienie dekompozycji spektralnej na macierze prostokątne,

Jeżeli Σ\Sigma jest symetryczna i dodatnio określona, to U=V=PU=V=P, więc C=Σ1/2=UD1/2VT, C = \Sigma^{1/2} = U D^{1/2} V^T, Funkcje fortranowe dla bibliotek LAPACK pozwalają na uwzględnienie faktu, że rozkładana macierz jest symetryczna i dodatnio określona, ale wrappery R już tego nie uwzględniają, więc ta transformacja jest mniej efektywna. Macierze UU, VV i DD dla macierzy Σ\Sigma można znaleźć z użyciem funkcji svd().

  • Dekompozycja Choleskiego, tzw. pierwiastek Choleskiego, Σ=QTQ, \Sigma = Q^T Q, gdzie QQ to macierz górna trójkątna. Macierz QQ dla macierzy Σ\Sigma można znaleźć z użyciem funkcji chol().

Który z tych algorytmów najlepiej wybrać? Szybszy i numerycznie bardziej stabilny! Porównajmy czasy działania różnych faktoryzacji.

  • Faktoryzacja (rozkład) Choleskiego.
n <- 200; d <- 100; N <- 10000
Sigma <- cov(matrix(rnorm(n*d),n,d))
X <- matrix(rnorm(n*d),n,d)
system.time(replicate(N, {Q <- chol(Sigma)
      X %*%  Q} ))
##    user  system elapsed 
##  20.557   2.025  23.770
  • Faktoryzacja SVD
system.time(replicate(N, {tmp <- svd(Sigma)
      X %*% (tmp$u %*% diag(sqrt(tmp$d)) %*% t(tmp$v))}))
##    user  system elapsed 
## 101.648   8.117 119.463
  • Rozkład spektralny / na wartości własne.
system.time(replicate(N, {tmp <- eigen(Sigma, symmetric=T)
      X %*% (tmp$vectors %*% diag(sqrt(tmp$values)) %*% t(tmp$vectors))
            }))
##    user  system elapsed 
##  66.521   5.227  78.655

Powyższe czasy zostały policzone dla pewnych przykładowych danych i mogą różnić się w zależności od zainstalowanych bibliotek. Ogólnie jednak najszybszą procedurą jest faktoryzacja Choleskiego i to ona najczęściej jest wykorzystywana do generowania zmiennych o wielowymiarowym rozkładzie normalnym.

W programie R nie musimy do generacji zmiennych z rozkładu normalnego wielowymiarowego ręcznie wykonywać dekompozycji macierzy Σ\Sigma, wystarczy użyć funkcji rmvnorm(mvtnorm). Pozwala ona wskazać metodę faktoryzacji argumentem method=c("eigen", "svd", "chol"). Domyślnie wykorzystywana jest dekompozycja na wartości spektralne (na wypadek gdyby macierz Σ\Sigma była niepełnego rzędu), ale jeżeli zależy nam na szybkości i mamy pewność że macierz Σ\Sigma jest dodatnio określona to lepiej użyć pierwiastka Choleskiego.