Statistische Informationsverarbeitung - Projekt

From StatWiki

Jump to: navigation, search

Contents


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





 Numbers of failures and times of observation of 10 pumps in a nuclear plant (Source: Gaver and O'Muircheartaigh 1987) 

     Pump 1   Pump 2   Pump 3   Pump 4   Pump 5   Pump 6   Pump 7   Pump 8   Pump 9   Pump 10   

  
  Failures   5.00   1.00   5.00   14.00   3.00   19.00   1.00   1.00   14.00   22.00  

  
  Time   94.32   15.72   62.88   125.76   5.24   31.44   1.05   1.05   2.10   10.48

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

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.



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

#


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

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



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.



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.



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

#
Iterations = 1:1000
Thinning interval = 1 
Number of chains = 9 
Sample size per chain = 1000 

1. Empirical mean and standard deviation for each variable,



  plus standard error of the mean:



           Mean      SD Naive SE Time-series SE



lambda[2] 0.1587 0.09526 0.001004       0.001059
lambda[8] 1.0164 0.63624 0.006707       0.007542

2. Quantiles for each variable:



            2.5%      5%    95%  97.5%



lambda[2] 0.03045 0.03996 0.3416 0.3951
lambda[8] 0.18610 0.25169 2.2315 2.6026










Retrieved from "http://www.eisber.net/StatWiki/index.php/Statistische_Informationsverarbeitung_-_Projekt"