@@ -0,0 +1,2103 @@
+import numpy as np%0Aimport pandas as pd%0A%0A# code adapted from pymc3.diagnostics:%0A# https://github.com/pymc-devs/pymc3/blob/master/pymc3/diagnostics.py%0A%0Adef _vhat_w(mcmc):%0A    # Calculate between-chain variance%0A    B = mcmc.n_samples * mcmc.data.mean(axis=0).var(axis=0, ddof=1)%0A%0A    # Calculate within-chain variance%0A    W = mcmc.data.var(axis=0, ddof=1).mean(axis=0)%0A%0A    # Estimate of marginal posterior variance%0A    Vhat = W * (mcmc.n_samples - 1) / mcmc.n_samples + B / mcmc.n_samples%0A%0A    return Vhat, W%0A%0Adef gelman_rubin(mcmc):%0A    %22%22%22%0A    mcmc (MCMCResults): Pre-sliced MCMC samples to compute diagnostics for.%0A    %22%22%22%0A%0A    if mcmc.n_chains %3C 2:%0A        raise ValueError(%0A            'Gelman-Rubin diagnostic requires multiple chains '%0A            'of the same length.')%0A%0A    # get Vhat and within-chain variance%0A    Vhat, W = _vhat_w(mcmc)%0A%0A    # compute and return Gelman-Rubin statistic%0A    Rhat = np.sqrt(Vhat / W)%0A    return pd.DataFrame(%7B'gelman_rubin':Rhat%7D, index=mcmc.levels)%0A%0Adef effective_n(mcmc):%0A    %22%22%22%0A    mcmc (MCMCResults): Pre-sliced MCMC samples to compute diagnostics for.%0A    %22%22%22%0A%0A    if mcmc.n_chains %3C 2:%0A        raise ValueError(%0A            'Calculation of effective sample size requires multiple chains '%0A            'of the same length.')%0A%0A    def get_neff(x, Vhat):%0A        n_samples = x.shape%5B0%5D%0A        n_chains = x.shape%5B1%5D%0A%0A        negative_autocorr = False%0A        t = 1%0A%0A        rho = np.ones(n_samples)%0A        # Iterate until the sum of consecutive estimates of autocorrelation is%0A        # negative%0A        while not negative_autocorr and (t %3C n_samples):%0A            variogram = np.mean((x%5Bt:, :%5D - x%5B:-t, :%5D)**2)%0A            rho%5Bt%5D = 1. - variogram / (2. * Vhat)%0A            negative_autocorr = sum(rho%5Bt - 1:t + 1%5D) %3C 0%0A            t += 1%0A%0A        if t %25 2:%0A            t -= 1%0A%0A        return min(n_chains * n_samples,%0A            int(n_chains * n_samples / (1. + 2 * rho%5B1:t-1%5D.sum())))%0A%0A    Vhat, W = _vhat_w(mcmc)%0A    n_eff = %5Bget_neff(mcmc.data%5B:, :, i%5D, x) for i, x in enumerate(Vhat)%5D%0A%0A    return pd.DataFrame(%7B'effective_n':n_eff%7D, index=mcmc.levels)%0A

Add diagnostics.py module