Statistische Informationsverarbeitung - Übung 6

From StatWiki
Jump to: navigation, search

11. A random walk with absorbing barriers

Sei X_n\, eine Irrfahrt auf X = \{0,1,\ldots,N\}\, bei der Sie zu jedem Zeitschritt mit Wahrscheinlichkeit p einen Schritt nach links (-1) und mit Wahrscheinlichkeit r ein Schritt nach rechts (+1) machen, oder aber mit Wahrscheinlichkeit q auf der Stelle treten.

(Es gilt p + q + r = 1.) Erreichen Sie allerdings einen der beiden "Ränder" (0 bzw. N) machen Sie mit Sicherheit einen Schritt zurück (also +1 bzw. -1).

(a)

Schreiben Sie die Übergangsmatrix \mathbb{P}\, von X_n\, auf.


 P = \begin{bmatrix}
0 & 1 & 0 & 0 & 0 & \ldots & 0 \\
p & q & r & 0 & 0 & \ldots & 0 \\
0 & p & q & r & 0 & \ldots & 0 \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
0 & \ldots & 0 & 0 & p & q & r \\
0 & \ldots & 0 & 0 & 0 & 1 & 0 
\end{bmatrix}

(b)

Berechne Sie für N = 3, p = r = \frac 1 4 \, und  q = \frac 1 2\, die stationäre Verteilung \vec \pi = (\pi_0, \pi_1, \pi_2, \pi_3) \,.

\pi = \pi P\,

\begin{align}
\pi_0 & = \frac 1 4 \pi_1 \\
\pi_3 & = \frac 1 4 \pi_2 \\
\pi_1 & = \pi_0 +  \frac 1 2 \pi_1 + \frac 1 4 \pi_2 \\
\pi_2 & = \frac 1 4 \pi_1 + \frac 1 2 \pi_2 + \pi_3 \\
1 & = \pi_0 + \pi_1 + \pi_2 + \pi_3 \\
\end{align}
 
\Rightarrow \vec \pi = \left(\frac{1}{10}, \frac{4}{10}, \frac{4}{10}, \frac{1}{10}\right)\,

(c)

Ist X_n\, ergodisch?

Ja, weil die Markovkette aus einer Klasse besteht und aperiodisch ist. (Zeichnung)

(d)

MCMC!

Generieren Sie \vec X \sim \vec \pi \, indem Sie sich auf die (oben definierte) Irrfahrt begeben.

Starten Sie in einem x_0 \in \mathcal \, Ihrer Wahl starten und führen Sie dann n Zeitschritte durch.

Wiederholen Sie diesen Vorgang jeweils, um 1000 und 10000 X's, also Werte in \mathcal X = \{0,1,2,3\} \, zu erzeugen.

Probieren Sie jeweils n = 10 und n = 100 Zeitschritte. Vergleichen Sie Ihre Ergebnisse (die relativen Häufigkeiten der Werte in \mathcal X = \{0,1,2,3\} \,) mit dem in (b) errechneten \vec \pi \,.

REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/99c3d081c9a4bab5db65d224b394cb3956b79ece_%i.pdf'
> rpdfno<-0
> rhtml<-''
> rfiles<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/'
> source('/var/www/localhost/htdocs/StatWiki/Rfiles/R/@.R')
> rout<-'text'
> cat('<!--- Start of program --->\n')
<!--- Start of program --->
> library(distrEx)
Error in library(distrEx) : there is no package called ‘distrEx’
Execution halted
in

library(distrEx) distroptions("WarningSim" = FALSE)

p=matrix(c(0,1,0,0,1/4,1/2,1/4,0,0,1/4,1/2,1/4,0,0,1,0),nrow=4,byrow=T)

finv=matrix(unlist(mapply(rep, 1:4, t(p*100))), byrow=T, ncol=100) p.dd = sapply(1:4,function(x) DiscreteDistribution(supp = c(1:4), prob = p[x,]) )

genpath = function(x0,n) { # I need 1 <= x <= 100 for the indicies ru = 1 + floor(runif(n,max=100)) # path = rep(1,n) # path[1] = x0 tab = rep(0,4) tab[x0] = tab[x0] + 1 for(j in 2:n) { # path[j] = finv[path[j-1], ru[j]] x0 = finv[x0, ru[j]] # x0 = r(p.dd[[x0]])(1) tab[x0] = tab[x0] + 1 } # path tab }

plot.res = function(res,main) { # prop.table(rowSums(res)) d=data.frame(count=c(apply(res, 2, prop.table)),state=rep(1:4,1000)) boxplot(count~state, data=d, names=c("z 0","z 1","z 2","z 3"),main=main,ylim=c(0,1),

               ylab="Relative Frequency in State")

for(i in 1:4) segments(i-.5,statdist[i], i+.5, statdist[i], col=2) }

  1. stationary distribution

statdist=c(1/10,2/5,2/5,1/10) statdist.inv = unlist(mapply(rep, 1:4, statdist*10))

statdist.dd=DiscreteDistribution(supp = c(1:4), prob = statdist)

Die roten Balken sind jeweils die Werte der stationären Verteilung

n=10, x_0 = 1\, mit 1000 X's und 10000 X's sowie n=100, x_0 = 1\, mit 1000 X's und 10000 X's

REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/ad4a3c66163f1666bbefc5840a21c26caf117b5b_%i.pdf'
> rpdfno<-0
> rhtml<-''
> rfiles<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/'
> source('/var/www/localhost/htdocs/StatWiki/Rfiles/R/@.R')
> rout<-'display'
> cat('<!--- Start of program --->\n')
<!--- Start of program --->
> pdf(rpdf,width=12,height=12)
> par(mfrow=c(2,2))
> plot.res(replicate(1000,  genpath(1,10)), "n=10 und 1000 X's")
Error: could not find function "plot.res"
Execution halted
in
 

pdf(rpdf,width=12,height=12) par(mfrow=c(2,2)) plot.res(replicate(1000, genpath(1,10)), "n=10 und 1000 X's") plot.res(replicate(10000, genpath(1,10)), "n=10 und 10000 X's") plot.res(replicate(1000, genpath(1,100)), "n=100 und 1000 X's") plot.res(replicate(10000, genpath(1,100)), "n=100 und 10000 X's")

(e)

Was passiert, wenn Sie Ihre Irrfahrt mit X_0 \sim \vec \pi \, starten und einen Zeitschritt (n = 1) durchführen?

Simulieren Sie so 1000 (oder mehr) Werte. Wie sieht die Verteilung der generierten Zahlen aus? War dies zu erwarten?

REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/f5e250c2b511bfff01dbe32a31d4d4bdac565405_%i.pdf'
> rpdfno<-0
> rhtml<-''
> rfiles<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/'
> source('/var/www/localhost/htdocs/StatWiki/Rfiles/R/@.R')
> rout<-'display'
> cat('<!--- Start of program --->\n')
<!--- Start of program --->
> library(distrEx)
Error in library(distrEx) : there is no package called ‘distrEx’
Execution halted
in
 

library(distrEx)

  1. prop.table(rowSums(replicate(1000, genpath(statdist.inv[1+runif(1,max=10)],100) )))

replicate(1000, genpath(r(statdist.dd)(1),100))

pdf(rpdf) plot.res(replicate(1000, genpath(r(statdist.dd)(1),100)), "n=100 und 1000 X's")

Die Verteilung ist sehr ähnlich zu \vec \pi\,.

Ja, weil sich \vec \pi \, mit einen neuem Schritt nicht verändert \left( = \vec \pi \mathbb P = \vec \pi \right) \,.

Metropolis-Hastings

Berechnen Sie die Übergangswahrscheinlichkeiten \alpha(x, y)\, für den Metropolis-Hastings Algorithmus mit target density \pi(x)\, und

ACHTUNG Buchstaben sind hier vertauscht!!!
Laut Definition ist  \alpha(x_{t},y_{t+1}) = \begin{cases} 
\min \left(\cfrac{\pi(y_{t+1})}{\pi(x_{t})} \cfrac{q(x_{t},y_{t+1})}{q(x_{t+1},y_{t})},1 \right) \qquad \mbox{falls } x_t \ne 0 \\ 
\mbox{sonst } 0
\end{cases} 

mit q(x,y) = q(x|y)\ ,

(a) einer symmetrischen candidate density q(x|y) \equiv q(y|x) \,

\begin{align}
 \alpha(x_{t},y_{t+1}) 
& = \min \left(\cfrac{\pi(y_{t+1})}{\pi(x_{t})} \cfrac{q(x_{t},y_{t+1})}{q(x_{t+1},y_{t})},1 \right) \\
& = \min \left(\cfrac{\pi(y_{t+1})}{\pi(x_{t})} \cfrac{q(x_{t},y_{t+1})}{q(x_{t},y_{t+1})},1 \right) \\
& = \min \left(\cfrac{\pi(y_{t+1})}{\pi(x_{t})},1 \right) 
\end{align}\, 

für x_t \ne 0\, und 0 sonst

(b) einer unabängigen candidate density q(y|x) \equiv q(y) \,

\begin{align}
 \alpha(x_{t},y_{t+1}) 
& = \min \left(\cfrac{\pi(y_{t+1})}{\pi(x_{t})} \cfrac{q(x_{t},y_{t+1})}{q(x_{t+1},y_{t})},1 \right) \\
& = \min \left(\cfrac{\pi(y_{t+1})}{\pi(x_{t})} \cfrac{q(x_{t})}{q(y_{t+1})},1 \right) 
\end{align}\, 

für x_t \ne 0\, und 0 sonst

Bemerkung (a) führt zum sog. Metropolis Algorithmus, (b) zum independence sampler.