Calculates Freeman-Tukey statistics (Freeman and Tukey 1950) as
    a measure of discrepancy between observed and r replicates of simulated data. This
    is a convenient method for assessing goodness-of-fit (see Brooks et al. 2000).

    D(x|\theta) = \sum_j (\sqrt{x_j} - \sqrt{e_j})^2

    :Parameters:
      observed : Iterable of observed values (size=(n,))
      simulated : Iterable of simulated values (size=(r,n))
      expected : Iterable of expected values (size=(r,) or (r,n))

    :Returns:
      D_obs : Discrepancy of observed values
      D_sim : Discrepancy of simulated values

def discrepancy(observed, simulated, expected):
    """Calculates Freeman-Tukey statistics (Freeman and Tukey 1950) as
    a measure of discrepancy between observed and r replicates of simulated data. This
    is a convenient method for assessing goodness-of-fit (see Brooks et al. 2000).

    D(x|\theta) = \sum_j (\sqrt{x_j} - \sqrt{e_j})^2

    :Parameters:
      observed : Iterable of observed values (size=(n,))
      simulated : Iterable of simulated values (size=(r,n))
      expected : Iterable of expected values (size=(r,) or (r,n))

    :Returns:
      D_obs : Discrepancy of observed values
      D_sim : Discrepancy of simulated values

    """
    try:
        simulated = simulated.astype(float)
    except AttributeError:
        simulated = simulated.trace().astype(float)
    try:
        expected = expected.astype(float)
    except AttributeError:
        expected = expected.trace().astype(float)
    # Ensure expected values are rxn
    expected = np.resize(expected, simulated.shape)

    D_obs = np.sum([(np.sqrt(observed) - np.sqrt(
        e)) ** 2 for e in expected], 1)
    D_sim = np.sum(
        [(np.sqrt(s) - np.sqrt(e)) ** 2 for s,
         e in zip(simulated,
                  expected)],
        1)

    # Print p-value
    count = sum(s > o for o, s in zip(D_obs, D_sim))
    print_('Bayesian p-value: p=%.3f' % (1. * count / len(D_obs)))

    return D_obs, D_sim