for ( i in 1:n ) { y.current <- y[ -i ] n.current <- length( y.current ) s.current <- sum( y.current ) alpha.star <- alpha + s.current beta.star <- beta + n.current predictive.mean.current <- alpha.star / beta.star predictive.SD.current <- sqrt( ( alpha.star / beta.star ) * ( 1 + 1 / beta.star ) ) z[ i ] <- ( y[ i ] - predictive.mean.current ) / predictive.SD.current }