Model validation method, following Cook et al. (Journal of Computational and
    Graphical Statistics, 2006, DOI: 10.1198/106186006X136976).

    Generates posterior samples based on 'true' parameter values and data simulated
    from the priors. The quantiles of the parameter values are calculated, based on
    the samples. If the model is valid, the quantiles should be uniformly distributed
    over [0,1].

    Since this relies on the generation of simulated data, all data stochastics
    must have a valid random() method for validation to proceed.

    Parameters
    ----------
    sampler : Sampler
      An MCMC sampler object.
    replicates (optional) : int
      The number of validation replicates (i.e. number of quantiles to be simulated).
      Defaults to 100.
    iterations (optional) : int
      The number of MCMC iterations to be run per replicate. Defaults to 2000.
    burn (optional) : int
      The number of burn-in iterations to be run per replicate. Defaults to 1000.
    thin (optional) : int
      The thinning factor to be applied to posterior sample. Defaults to 1 (no thinning)
    deterministic (optional) : bool
      Flag for inclusion of deterministic nodes in validation procedure. Defaults
      to False.
    db (optional) : string
      The database backend to use for the validation runs. Defaults to 'ram'.
    plot (optional) : bool
      Flag for validation plots. Defaults to True.

    Returns
    -------
    stats : dict
      Return a dictionary containing tuples with the chi-square statistic and
      associated p-value for each data stochastic.

    Notes
    -----
    This function requires SciPy.

def validate(sampler, replicates=20, iterations=10000, burn=5000,
             thin=1, deterministic=False, db='ram', plot=True, verbose=0):
    """
    Model validation method, following Cook et al. (Journal of Computational and
    Graphical Statistics, 2006, DOI: 10.1198/106186006X136976).

    Generates posterior samples based on 'true' parameter values and data simulated
    from the priors. The quantiles of the parameter values are calculated, based on
    the samples. If the model is valid, the quantiles should be uniformly distributed
    over [0,1].

    Since this relies on the generation of simulated data, all data stochastics
    must have a valid random() method for validation to proceed.

    Parameters
    ----------
    sampler : Sampler
      An MCMC sampler object.
    replicates (optional) : int
      The number of validation replicates (i.e. number of quantiles to be simulated).
      Defaults to 100.
    iterations (optional) : int
      The number of MCMC iterations to be run per replicate. Defaults to 2000.
    burn (optional) : int
      The number of burn-in iterations to be run per replicate. Defaults to 1000.
    thin (optional) : int
      The thinning factor to be applied to posterior sample. Defaults to 1 (no thinning)
    deterministic (optional) : bool
      Flag for inclusion of deterministic nodes in validation procedure. Defaults
      to False.
    db (optional) : string
      The database backend to use for the validation runs. Defaults to 'ram'.
    plot (optional) : bool
      Flag for validation plots. Defaults to True.

    Returns
    -------
    stats : dict
      Return a dictionary containing tuples with the chi-square statistic and
      associated p-value for each data stochastic.

    Notes
    -----
    This function requires SciPy.
    """
    import scipy as sp
    # Set verbosity for models to zero
    sampler.verbose = 0

    # Specify parameters to be evaluated
    parameters = sampler.stochastics
    if deterministic:
        # Add deterministics to the mix, if requested
        parameters = parameters | sampler.deterministics

    # Assign database backend
    original_backend = sampler.db.__name__
    sampler._assign_database_backend(db)

    # Empty lists for quantiles
    quantiles = {}

    if verbose:
        print_("\nExecuting Cook et al. (2006) validation procedure ...\n")

    # Loop over replicates
    for i in range(replicates):

        # Sample from priors
        for p in sampler.stochastics:
            if not p.extended_parents:
                p.random()

        # Sample "true" data values
        for o in sampler.observed_stochastics:
            # Generate simuated data for data stochastic
            o.set_value(o.random(), force=True)
            if verbose:
                print_("Data for %s is %s" % (o.__name__, o.value))

        param_values = {}
        # Record data-generating parameter values
        for s in parameters:
            param_values[s] = s.value

        try:
            # Fit models given parameter values
            sampler.sample(iterations, burn=burn, thin=thin)

            for s in param_values:

                if not i:
                    # Initialize dict
                    quantiles[s.__name__] = []
                trace = s.trace()
                q = sum(trace < param_values[s], 0) / float(len(trace))
                quantiles[s.__name__].append(open01(q))

            # Replace data values
            for o in sampler.observed_stochastics:
                o.revert()

        finally:
            # Replace data values
            for o in sampler.observed_stochastics:
                o.revert()

            # Replace backend
            sampler._assign_database_backend(original_backend)

        if not i % 10 and i and verbose:
            print_("\tCompleted validation replicate", i)

    # Replace backend
    sampler._assign_database_backend(original_backend)

    stats = {}
    # Calculate chi-square statistics
    for param in quantiles:
        q = quantiles[param]
        # Calculate chi-square statistics
        X2 = sum(sp.special.ndtri(q) ** 2)
        # Calculate p-value
        p = sp.special.chdtrc(replicates, X2)

        stats[param] = (X2, p)

    if plot:
        # Convert p-values to z-scores
        p = copy(stats)
        for i in p:
            p[i] = p[i][1]
        Matplot.zplot(p, verbose=verbose)

    return stats