Python(https://github.com/pymc-devs/pymc/blob/2.3/pymc/diagnostics.py)(核心部分):
<br />
for i in range(intervals):<br />
# Size of remaining array<br />
n = len(x_trunc)<br />
# Calculate slices<br />
first_slice = x_trunc[:int(first * n)]<br />
last_slice = x_trunc[int(last * n):]<br />
z = (first_slice.mean() - last_slice.mean())<br />
z /= np.sqrt(first_slice.var() / len(first_slice) + last_slice.var() / len(last_slice))<br />
zscores[i] = len(x) - len(x_trunc), z<br />
x_trunc = x_trunc[int(first * n):]<br />
return zscores<br />
</p>
R的coda包中的geweke.diag函数(完整):
<br />
> geweke.diag<br />
function (x, frac1 = 0.1, frac2 = 0.5)<br />
{<br />
if (frac1 < 0 || frac1 > 1) {<br />
stop("frac1 invalid")<br />
}<br />
if (frac2 < 0 || frac2 > 1) {<br />
stop("frac2 invalid")<br />
}<br />
if (frac1 + frac2 > 1) {<br />
stop("start and end sequences are overlapping")<br />
}<br />
if (is.mcmc.list(x)) {<br />
return(lapply(x, geweke.diag, frac1, frac2))<br />
}<br />
x <- as.mcmc(x)<br />
xstart <- c(start(x), floor(end(x) - frac2 * (end(x) - start(x))))<br />
xend <- c(ceiling(start(x) + frac1 * (end(x) - start(x))), end(x))<br />
y.variance <- y.mean <- vector("list", 2)<br />
for (i in 1:2) {<br />
y <- window(x, start = xstart[i], end = xend[i])<br />
y.mean[[i]] <- apply(as.matrix(y), 2, mean)<br />
y.variance[[i]] <- spectrum0.ar(y)$spec/niter(y)<br />
}<br />
z <- (y.mean[[1]] - y.mean[[2]])/sqrt(y.variance[[1]] + y.variance[[2]])<br />
out <- list(z = z, frac = c(frac1, frac2))<br />
class(out) <- "geweke.diag"<br />
return(out)<br />
}<br />
</p>