code: def discrepancy(observed, simulated, expected):
      
      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
docs: 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