Statistische Informationsverarbeitung - Projekt

From StatWiki
Jump to: navigation, search

Nuclear pump failures

Gaver and O'Muircheartaigh (1987) introduced a model that is frequently used (or even overused) in the Gibbs sampling literature to illustrate various properties (see, for instance, Gelfand and Smith 1990, Tanner 1996, or Guihenneuc-Jouyaux and Robert 1998).

Engine.php: unknown value "raw" for attribute "output"
in

pumps=matrix(c(5,1,5,14,3,19,1,1,14,22,94.32,15.72,62.88,125.76,5.24,31.44,1.05,1.05,2.10,10.48),byrow=T,nrow=2) rownames(pumps) = rownames=c("Failures", "Time") colnames(pumps) = paste("Pump",1:10)

library(xtable) print(xtable(pumps,caption="Numbers of failures and times of observation of 10 pumps in a nuclear plant (Source: Gaver and O'Muircheartaigh 1987)"), type="html")

This model describes multiple failures of pumps in a nuclear plant, with the data given in the table above.

The modeling is based on the assumption that the failures of the i-th pump follow a Poisson process with parameters \lambda_i (1 \le i \le 10) \,.

For an observed time t_i\,, the number of failures p_i\, is thus a Poisson P(\lambda_i, t_i)\, random variable. The associated prior distributions are (1 \le i \le 10) \,

\lambda_i \sim \mbox{iid Gamma}(\alpha, \beta), \qquad \beta \sim \mbox{Gamma}(\gamma, \delta) \,

with \alpha = 1.8, \gamma = 0.01, \delta = 1\,

Beispiel (a)

Berechnen Sie den Posterior \pi(\lambda,\beta|X)\,, wobei \lambda = (\lambda_1, \ldots, \lambda_n)\, und n=10\, ist. \alpha, \gamma, \delta \, sind fix.

 \mbox{Posterior} \propto \mbox{Likelihood} \times \mbox{Prior} 
Laut Angabe:

 \begin{align}
f(\lambda_i|\beta)
& = \Gamma(\alpha)^{-1} \beta^\alpha\lambda_i^{\alpha-1}e^{-\beta\lambda_i}\\
\end{align}
 
f(\beta) = \Gamma(\gamma)^{-1} \delta^{\gamma}\beta^{\gamma-1}e^{-\delta\beta} \,
Satz der bedingten Wahrscheinlichkeit und weil \lambda_i\, i.i.d. sind

 \begin{align}
\mbox{Prior} 
& = f(\lambda_1,\ldots,\lambda_n,\beta)\\
& = f(\lambda_1, \ldots ,\lambda_n|\beta) \cdot f(\beta)\\
& = f(\lambda_1|\beta) \cdots f(\lambda_n|\beta) \cdot f(\beta)\\
& = \prod_{i=1}^nf(\lambda_i|\beta) \cdot f(\beta)
\end{align} 

 \begin{align}
f(\lambda_1,\ldots,\lambda_n|\beta)
& = \prod_{i=1}^n f(\lambda_i|\beta) \\
& = \prod_{i=1}^n \Gamma(\alpha)^{-1}\beta^\alpha\lambda_i^{\alpha-1}e^{-\beta\lambda_i}\\
& = \Gamma(\alpha)^{-n} \beta^{n\cdot\alpha} \prod_{i=1}^n \lambda_i^{\alpha-1} e^{-\beta\lambda_i} \\
\end{align} 
 \begin{align}
\mbox{Prior} 
& = f(\beta) \prod_{i=1}^nf(\lambda_i|\beta) \\
& = 
\left[ \Gamma(\gamma)^{-1} \delta^{\gamma}\beta^{\gamma-1}e^{-\delta\beta} \right]
\left[ \Gamma(\alpha)^{-n} \beta^{n\cdot\alpha} \prod_{i=1}^n\lambda_i^{\alpha-1} e^{-\beta\lambda_i} \right] \\
& = \Gamma(\alpha)^{-n} \Gamma(\gamma)^{-1} \beta^{n\cdot\alpha + (\gamma-1)} \delta^{\gamma} e^{- \delta\beta}
\prod_{i=1}^n\lambda_i^{\alpha-1} e^{-\beta\lambda_i} \\
& \propto \beta^{n\cdot\alpha + (\gamma-1)} e^{- \delta\beta}
  \prod_{i=1}^n\lambda_i^{\alpha-1} e^{-\beta\lambda_i} \\
\end{align}

 \begin{align}
\mbox{Likelihood} 
& = \prod_{i=1}^n\frac{(\lambda_it_i)^{p_i}}{p_i!}e^{-\lambda_i t_i} \\
& \propto \prod_{i=1}^n \lambda_i^{p_i} e^{-\lambda_i t_i} 
\end{align}  

Hier werden die Daten eingesetzt, 
d.h. p_i\, sind die Failures und t_i\, die jeweiligen Zeiten. 
Die sind Konstant und können daher im nächsten Schritt vernachlässigt werden.

 \begin{align}
\mbox{Posterior} 
& \propto \mbox{Likelihood} \times \mbox{Prior} \\
& \propto 
  \left[ \prod_{i=1}^n \lambda_i^{p_i} e^{-\lambda_i t_i} \right]
  \left[ \beta^{n\cdot\alpha + (\gamma-1)} e^{- \delta\beta} \prod_{i=1}^n\lambda_i^{\alpha-1} e^{-\beta\lambda_i} \right]\\
& \propto  \beta^{n\alpha+\gamma-1}e^{-\beta\delta}\prod_{i=1}^n\lambda_i^{\alpha+p_i-1}e^{-\lambda_i(\beta+t_i)}
\end{align}

Beispiel (b)

Berechnen Sie die bedingten Verteilung für \lambda_i|\lambda_{i-1},\beta,X\, und i = 1, \ldots, n \,.

Dichte einer Gamma-Verteilung

\mbox{Gamma}(x)
= \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha - 1} e^{-x \cdot \beta} 
\propto x^{\alpha - 1} e^{-x \cdot \beta}
\,

 \begin{matrix}
\mbox{Posterior} \propto 
 & \underbrace{\lambda_1^{\alpha+p_1-1}} 
 & \underbrace{e^{-\lambda_1(\beta+t_1)}}
 & \underbrace{\beta^{n\alpha+\gamma-1}e^{-\beta\delta}  \prod_{i=2}^n\lambda_i^{\alpha+p_i-1}e^{-\lambda_i(\beta+t_i)}} \\
& x^{\alpha-1} & e^{-x \cdot \beta} & \mbox{some constant}
\end{matrix}

Alles, was nach dem Produktzeichen steht, ist ja nur eine Konstante.
 
Der Zusammenhang zur Dichte der Gammaverteilung ist nun durch x = \lambda_1, \alpha=\alpha + p_1, \beta=\beta+t_1\, gegeben, wobei

links immer die Variablen der Dichte und rechts die Variablen aus der Posterior stehen.
Daraus kann man ablesen (einfach die Funktion so umschreiben wie oben), dass die bedingte Dichte von \lambda_i\, ist: 

 \lambda_i \sim Ga(\alpha+p_i,\beta+t_i) \,

Beispiel (c)

Berechnen Sie die bedingten Verteilung für \beta|\lambda_{i},X\, und i = 1, \ldots, n \,.

 \mbox{Posterior} \propto \beta^{n \cdot \alpha+\gamma-1} e^{-\beta \left(\delta + \sum_{i=1}^n \lambda_i\right) } \prod_{i=1}^n\lambda_i^{\alpha+p_i-1}e^{-t_i} 

 \beta \sim Ga \left(n \cdot \alpha+\gamma,\delta+\sum_{i=1}^n\lambda_i \right) 

Beispiel (d)

Programmieren Sie einen Gibbs Sampler für die oben angegeben Werte von  \alpha, \gamma, \delta \,

Engine.php: unknown value "none" for attribute "output"
in

library(coda) alpha = 1.8 gamma = 0.01 delta = 1

r.lambda.i = function(beta, i) rgamma(1, alpha + pumps[1,i], beta + pumps[2,i]) r.beta = function(lambda) rgamma(1, length(lambda) * alpha + delta, delta + sum(lambda))

gibbs = function(beta,n=1000) {

d = matrix(0, nrow=n, ncol=11) colnames(d) = c("beta",sprintf("lambda[%d]", 1:10))

for(i in 1:n) { for(j in 1:10) d[i,j+1] = r.lambda.i(beta,j) d[i,1] = beta = r.beta(d[i,-1]) }


       mcmc(d)

}

mc = mcmc.list(sapply(0:8, gibbs, simplify=F))

Beispiel (e)

Betracten Sie die fortlaufenden Mittel und den Gelman-Rubion Plot von  \lambda_2, \lambda_8, \beta \,, um eine geeignete Anzahl von Iterationen festzulegen.

REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/398f5ece6898ea22dd3c4494c1c52955e2520283_%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,height=12, width=12)
> library(coda)
Loading required package: lattice
> gelman.plot(mc[,c(1,3,9)])
Error in as.mcmc.list(x) : object 'mc' not found
Calls: gelman.plot -> as.mcmc.list
Execution halted
in

pdf(rpdf,height=12, width=12) library(coda) gelman.plot(mc[,c(1,3,9)])

REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/a08a88ba677232dce8823d498cdc93f55dea17e2_%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,height=12,width=12)
> library(coda)
Loading required package: lattice
>
> # matrix an fortlaufenden mittel
> plot.mean = function(var.index) {
+ 	fm = matrix(unlist(lapply(mc[,var.index], function(x) cumsum(x)/(1:length(x)))), ncol=9)
+
+ 	var.name = varnames(mc)[var.index]
+ 	matplot(fm, type="l",xlab="Iterationen",ylab=paste("Fortlaufendes Mittel ",var.name))
+ 	title(paste("Fortlaufendes Mittel von",var.name))
+ }
>
> par(mfrow=c(2,2))
> plot.mean(1)
Error in lapply(mc[, var.index], function(x) cumsum(x)/(1:length(x))) :
  object 'mc' not found
Calls: plot.mean -> matrix -> unlist -> lapply
Execution halted
in

pdf(rpdf,height=12,width=12) library(coda)

  1. matrix an fortlaufenden mittel

plot.mean = function(var.index) { fm = matrix(unlist(lapply(mc[,var.index], function(x) cumsum(x)/(1:length(x)))), ncol=9)

var.name = varnames(mc)[var.index] matplot(fm, type="l",xlab="Iterationen",ylab=paste("Fortlaufendes Mittel ",var.name)) title(paste("Fortlaufendes Mittel von",var.name)) }

par(mfrow=c(2,2)) plot.mean(1) plot.mean(3) plot.mean(9)

Aus den oben angführten Plots ist zu erkennen, dass ab 400 Iterationen ein recht gutes Ergebnis erzielt werden kann.

Beispiel (f)

Plotten Sie 95% Konfidenzbänder für  \lambda_2 \, und  \lambda_8 \,, zusammen mit den entsprechenden fortlaufenden Mitteln.

REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/6e46e57471a6c6895ae633efde0caad77dac4d6e_%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)
> library(coda)
Loading required package: lattice
> plot.glob.mean = function(var.index) {
+ 	r = matrix(unlist(mc[,var.index]), nrow=length(mc))
+ 	r = apply(r, 2, mean)
+ 	fm = cumsum(r)/(1:length(r))
+ 	co = t(sapply(1:length(r), function(t) quantile(r[1:t], c(0.025,0.975))))
+
+ 	# mean
+ 	var.name = varnames(mc)[var.index]
+ 	plot(fm, type="l",xlab="Iterationen",ylab=paste("Gesamt fortlaufendes Mittel ",var.name), ylim=range(co))
+ 	title(paste("Gesamt fortlaufendes Mittel von",var.name))
+
+ 	# confidence interval
+ 	matlines(co, type="l", col=2, lty=2)
+ }
>
> par(mfrow=c(1,2))
> plot.glob.mean(3)
Error in unlist(mc[, var.index]) : object 'mc' not found
Calls: plot.glob.mean -> matrix -> unlist
Execution halted
in

pdf(rpdf,width=12) library(coda) plot.glob.mean = function(var.index) { r = matrix(unlist(mc[,var.index]), nrow=length(mc)) r = apply(r, 2, mean) fm = cumsum(r)/(1:length(r)) co = t(sapply(1:length(r), function(t) quantile(r[1:t], c(0.025,0.975))))

# mean var.name = varnames(mc)[var.index] plot(fm, type="l",xlab="Iterationen",ylab=paste("Gesamt fortlaufendes Mittel ",var.name), ylim=range(co)) title(paste("Gesamt fortlaufendes Mittel von",var.name))

# confidence interval matlines(co, type="l", col=2, lty=2) }

par(mfrow=c(1,2)) plot.glob.mean(3) plot.glob.mean(9)

Alternativ auch direkt mit Funktionen aus dem coda Package

Achtung: hier wird nur die erste Kette verwendet und nicht wie oben, über alle Ketten gemittelt.

Hier sieht man auch recht gut, dass durch das mitteln über alle Ketten die Konfidenzbänder um einiges enger werden.

REngine.php: > rpdf<-'/var/www/localhost/htdocs/StatWiki/Rfiles/R/63a3672d2f60520312093894b0c029e609966d80_%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)
> library(coda)
Loading required package: lattice
>
> cumuplot(mc[[1]][,c(3)],ylab="lambda[2]",col=c(2,1))
Error in inherits(x, "mcmc") : object 'mc' not found
Calls: cumuplot -> par -> set.mfrow -> nvar -> is.mcmc -> inherits
Execution halted
in

pdf(rpdf) library(coda)

cumuplot(mc[[1]][,c(3)],ylab="lambda[2]",col=c(2,1)) cumuplot(mc[[1]][,c(9)],ylab="lambda[8]",col=c(2,1))

Beispiel (g)

Geben Sie Ihren MCMC-Schätzer für  \lambda_2 \, und  \lambda_8 \, (also die Störungsrate für den Poisson Prozess der 2. und 8. Pumpe) jeweils mit einem 90% und 95% Konfidenzintervall, an.

Engine.php: unknown value "raw" for attribute "output"
in

library(coda) summary(mc[,c(3,9)],quantiles=c(0.025,0.05,.95,.975))