1. Εάν \(X\sim G(α,λ)\), τότε η ροπογεννήτριά της είναι: \(Μ_Χ(t) = E[e^{tX}] = \int_{0}^{\infty} e^{tx}\frac{λ^α}{Γ(α)}x^{α-1}e^{-λx} \,dx = \frac{λ^α}{Γ(α)}\int_{0}^{\infty} e^{(t-λ)x}x^{α-1} \,dx\). Εάν θέσουμε \(u = x(λ-t)\), το ολοκλήρωμα γίνεται: \(\frac{λ^α}{Γ(α)}\int_{0}^{\infty} e^{-u}(\frac{u}{λ-t})^{α-1} \,\frac{dx}{λ-t} = \frac{λ^α}{Γ(α)}\frac{1}{(λ-t)^α}\int_{0}^{\infty} e^{-u}u^{α-1} \,dx = (\frac{λ}{λ-t})^a,\)όπου \(t<λ\). Επίσης, εάν \(Y\sim Exp(λ)\), η ροπογεννήτριά της δίνεται από: \(M_Y(t) = \int_{0}^{\infty}e^{ty}λe^{-λt}\,dy = λ\int_{0}^{\infty}e^{-(λ-t)y}\,dy = \frac{λ}{λ-t}\int_{0}^{\infty}(λ-t)e^{-(λ-t)y}\,dy = \frac{λ}{λ-t},\) όπου \(t<λ\). Παρατηρούμε, λοιπόν, ότι αν \(S = \displaystyle\sum_{i=1}^{α}Y_i,\) με \(Y_i, Y_j\sim Exp(λ)\) ανεξάρτητες \(\forall i\neq j\), τότε κάθε μία από αυτές έχει ροπογεννήτρια \(M_{Y_i}(t)\) και \(M_S(t) = (M_Y(t))^α = Μ_{Y_1}\cdots M_{Y_α} = Μ_{Y_1+\dots+Y_α}\), αφού \(Y_1,\dots,Y_α\) ανεξάρτητες και ισόνομες. Δηλαδή, \(S\sim G(α,λ)\), και η ροπογεννήτρια μιας τ.μ. χαρακτηρίζει μονοσήμαντα την κατανομή της, επομένως πράγματι η \(X = \displaystyle\sum_{i=1}^αY_i\).
  1. Αρχικά, αν \(f(x)\) η σ.π.π. της \(G(α,λ_0)\), θεωρούμε υποψήφια πυκνότητα \(g(x)\) την σ.π.π. της \(Exp(λ)\) και υπολογίζουμε τον λόγο \(\frac{f(x)}{g(x)} = \frac{λ_0^α}{λΓ(α)}x^{α-1}e^{-(λ_0 - λ)x}\). Ξέρουμε ότι \((0,\infty)=S_f \subset S_g = [0,\infty)\) και θέλουμε, επίσης: \(\displaystyle\sup_{x\geq0}\frac{f(x)}{g(x)} = \displaystyle \sup_{x\geq0} x^{α-1}e^{-(λ_0-λ)x}< +\infty \Leftrightarrow λ\leqλ_0\), αφού \(0<α \in \mathbb{Z}\). Θέλουμε, τώρα, να βρούμε το \(Μ= \displaystyle\sup_{x\geq0}\frac{f(x)}{g(x)}\), οπότε αν \(h(x) = x^{α-1}e^{-(λ_0-λ)x}\) παραγωγίζουμε για να βρούμε το sup: \(h'(x) = 0 \Leftrightarrow x^{α-2}e^{-(λ_0-λ)x}[(α-1)-(λ_0-λ)x]=0\Leftrightarrow x=0 \;(α>2) \;ή\;x=\frac{α-1}{λ_0-λ}\geq0\). Tο \(x=0\) δεν μας δίνει supremum. Ας εξετάσουμε το \(x=\frac{α-1}{λ_0-λ}\), \(λ\neqλ_0\): παρατηρούμε ότι \(h''(\frac{α-1}{λ-λ_0})=-(λ_0-λ)(\frac{α-1}{λ_0-λ})^{α-2}e^{-(α-1)}<0 ,α\neq1\;λ\neqλ_0\), άρα είναι το supremum της \(h(x)\). Αντικαθιστούμε: \(Μ = Μ(λ) = \frac{λ_0^α}{λΓ(α)}(\frac{α-1}{λ_0-λ})^{α-1}e^{-(α-1)},\;α\neq1, λ_0\neqλ\) και βρίσκουμε το \(\displaystyle\sup_{0<λ<λ_0}λ(λ_0-λ)^{α-1}=\frac{λ_0}{α}>0\), διότι είναι το μοναδικό στάσιμο σημείο της \(K(λ) = λ(λ_0-λ)^{α-1}, λ<λ_0, α\neq1\) και \(K''(λ)= -4\frac{λ_0^2(α-1)^{α-1}}{α} <0, για\; α\neq1\).Τελικά, εάν \(λ\neqλ_0, \displaystyle\inf_{0<λ<λ_0}Μ(λ)=\frac{αλ_0^{α-1}}{Γ(α)}(\frac{α}{λ_0})^{α-1}e^{-(α-1)}\). Ας εξετάσουμε και την περίπτωση \(λ=λ_0\): τότε \(\displaystyle\sup_{0\leq x\leq1}x^{α-1} = 1 \;για\; x=1\) και \(Μ(λ_0)=\frac{λ_0^{α-1}}{Γ(α)}\). Παρατηρούμε ότι \(\frac{Μ(λ)}{Μ(λ_0)}=\frac{α^α}{(λ_0e)^{α-1}}>1, για\; α\neq1\), άρα η μέγιστη πιθανότητα αποδοχής (δηλαδή το μικρότερο Μ) επιτυγχάνεται για \(λ=λ_0\). Επίσης, για \(α=1\) έχουμε \(\displaystyle\sup_{x\geq0}e^{-(λ_0-λ)x} = 1, αφού\; λ\leqλ_0, άρα\; x=0\). Επομένως, θέλουμε \(\displaystyle\inf_{0<λ\leqλ_0}Μ(λ)=\frac{λ_0^α}{Γ(α)}\displaystyle\inf_{0<λ\leqλ_0}\frac{1}{λ}=\frac{λ_0^{α-1}}{Γ(α)}\), άρα και πάλι καταλήγουμε ότι η μέγιστη πιθανότητα αποδοχής δίνεται όταν \(λ=λ_0\).
sumexp <-function(nsim){
  y1<-rexp(nsim,2)
  y2<-rexp(nsim,2)
  y3<-rexp(nsim,2)
  y<-y1+y2+y3
  y
}
y<-sumexp(1000)
hist(y,main="Gamma(3,2) sample from Exponential",freq=F,col="grey")
curve((dgamma(x, shape = 3, rate = 2)), from = 0, to = 6, lty= 1, lwd = 1, col = "magenta", add=T)

Για να εφαρμόσουμε τον αλγόριθμο Αποδοχής-Απόρριψης, υπολογίζουμε πρώτα το \(Μ=Μ(2)=\frac{2^{3-1}}{ Γ(3)}=\frac{4}{2!}=2\).

acc.rej<-function(Nsim){
  M<-2
  x<-NULL
  while (length(x)<Nsim){
   u<-runif(Nsim*M)
   x<-c(x,y[u*M*dexp(Nsim*M,2)<dgamma(y, shape = 3, rate = 2)])
}
x<-x[1:Nsim]
x
}
y<-acc.rej(1000)
hist(y,main="Gamma(3, 2) sample with Accept-Reject",freq=F,col="grey")
curve((dgamma(x, shape = 3, rate = 2)),  from = 0, to = 6, lty= 1, lwd = 1, col = "magenta", add=T)

library(microbenchmark)
microbenchmark(sumexp(1000), acc.rej(1000), times=100)
## Unit: microseconds
##           expr   min     lq     mean median      uq     max neval
##   sumexp(1000) 586.4  644.6  928.110  665.1  807.60 20434.8   100
##  acc.rej(1000) 975.0 1063.5 1192.361 1141.1 1287.15  1827.2   100

Παρατηρούμε, λοιπόν, ότι εδώ η χρήση της μεθόδου Αποδοχής-Απόρριψης παράγει κατά μέσο όρο πιο αργά το ζητούμενο δείγμα από την προσομοίωση μέσω Εκθετικής κατανομής!