Statistische Informationsverarbeitung - Übung 6

From StatWiki

Jump to: navigation, search

Contents

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 \,. 

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)
}

# 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



 
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?



 
library(distrEx)

# 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.







Retrieved from "http://www.eisber.net/StatWiki/index.php/Statistische_Informationsverarbeitung_-_%C3%9Cbung_6"