#NEWLINE# def powerdiscrepancy( observed , expected , lambd = 0.0 , axis = 0 , ddof = 0 ) : #NEWLINE# #INDENT# "Calculates power discrepancy, a class of goodness-of-fit tests\n    as a measure of discrepancy between observed and expected data.\n\n    This contains several goodness-of-fit tests as special cases, see the\n    describtion of lambd, the exponent of the power discrepancy. The pvalue\n    is based on the asymptotic chi-square distribution of the test statistic.\n\n    freeman_tukey:\n    D(x|\theta) = \\sum_j (\\sqrt{x_j} - \\sqrt{e_j})^2\n\n    Parameters\n    ----------\n    o : Iterable\n        Observed values\n    e : Iterable\n        Expected values\n    lambd : float or string\n        * float : exponent `a` for power discrepancy\n        * 'loglikeratio': a = 0\n        * 'freeman_tukey': a = -0.5\n        * 'pearson': a = 1   (standard chisquare test statistic)\n        * 'modified_loglikeratio': a = -1\n        * 'cressie_read': a = 2/3\n        * 'neyman' : a = -2 (Neyman-modified chisquare, reference from a book?)\n    axis : int\n        axis for observations of one series\n    ddof : int\n        degrees of freedom correction,\n\n    Returns\n    -------\n    D_obs : Discrepancy of observed values\n    pvalue : pvalue\n\n\n    References\n    ----------\n    Cressie, Noel  and Timothy R. C. Read, Multinomial Goodness-of-Fit Tests,\n        Journal of the Royal Statistical Society. Series B (Methodological),\n        Vol. 46, No. 3 (1984), pp. 440-464\n\n    Campbell B. Read: Freeman-Tukey chi-squared goodness-of-fit statistics,\n        Statistics & Probability Letters 18 (1993) 271-278\n\n    Nobuhiro Taneichi, Yuri Sekiya, Akio Suzukawa, Asymptotic Approximations\n        for the Distributions of the Multinomial Goodness-of-Fit Statistics\n        under Local Alternatives, Journal of Multivariate Analysis 81, 335?359 (2002)\n    Steele, M. 1,2, C. Hurst 3 and J. Chaseling, Simulated Power of Discrete\n        Goodness-of-Fit Tests for Likert Type Data\n\n    Examples\n    --------\n\n    >>> observed = np.array([ 2.,  4.,  2.,  1.,  1.])\n    >>> expected = np.array([ 0.2,  0.2,  0.2,  0.2,  0.2])\n\n    for checking correct dimension with multiple series\n\n    >>> powerdiscrepancy(np.column_stack((observed,observed)).T, 10*expected, lambd='freeman_tukey',axis=1)\n    (array([[ 2.745166,  2.745166]]), array([[ 0.6013346,  0.6013346]]))\n    >>> powerdiscrepancy(np.column_stack((observed,observed)).T, 10*expected,axis=1)\n    (array([[ 2.77258872,  2.77258872]]), array([[ 0.59657359,  0.59657359]]))\n    >>> powerdiscrepancy(np.column_stack((observed,observed)).T, 10*expected, lambd=0,axis=1)\n    (array([[ 2.77258872,  2.77258872]]), array([[ 0.59657359,  0.59657359]]))\n    >>> powerdiscrepancy(np.column_stack((observed,observed)).T, 10*expected, lambd=1,axis=1)\n    (array([[ 3.,  3.]]), array([[ 0.5578254,  0.5578254]]))\n    >>> powerdiscrepancy(np.column_stack((observed,observed)).T, 10*expected, lambd=2/3.0,axis=1)\n    (array([[ 2.89714546,  2.89714546]]), array([[ 0.57518277,  0.57518277]]))\n    >>> powerdiscrepancy(np.column_stack((observed,observed)).T, expected, lambd=2/3.0,axis=1)\n    (array([[ 2.89714546,  2.89714546]]), array([[ 0.57518277,  0.57518277]]))\n    >>> powerdiscrepancy(np.column_stack((observed,observed)), expected, lambd=2/3.0, axis=0)\n    (array([[ 2.89714546,  2.89714546]]), array([[ 0.57518277,  0.57518277]]))\n\n    each random variable can have different total count/sum\n\n    >>> powerdiscrepancy(np.column_stack((observed,2*observed)), expected, lambd=2/3.0, axis=0)\n    (array([[ 2.89714546,  5.79429093]]), array([[ 0.57518277,  0.21504648]]))\n    >>> powerdiscrepancy(np.column_stack((observed,2*observed)), expected, lambd=2/3.0, axis=0)\n    (array([[ 2.89714546,  5.79429093]]), array([[ 0.57518277,  0.21504648]]))\n    >>> powerdiscrepancy(np.column_stack((2*observed,2*observed)), expected, lambd=2/3.0, axis=0)\n    (array([[ 5.79429093,  5.79429093]]), array([[ 0.21504648,  0.21504648]]))\n    >>> powerdiscrepancy(np.column_stack((2*observed,2*observed)), 20*expected, lambd=2/3.0, axis=0)\n    (array([[ 5.79429093,  5.79429093]]), array([[ 0.21504648,  0.21504648]]))\n    >>> powerdiscrepancy(np.column_stack((observed,2*observed)), np.column_stack((10*expected,20*expected)), lambd=2/3.0, axis=0)\n    (array([[ 2.89714546,  5.79429093]]), array([[ 0.57518277,  0.21504648]]))\n    >>> powerdiscrepancy(np.column_stack((observed,2*observed)), np.column_stack((10*expected,20*expected)), lambd=-1, axis=0)\n    (array([[ 2.77258872,  5.54517744]]), array([[ 0.59657359,  0.2357868 ]]))\n\n\n    " #NEWLINE# o = np . array ( observed ) #NEWLINE# e = np . array ( expected ) #NEWLINE# if ( np . isfinite ( lambd ) == True ) : #NEWLINE# #INDENT# a = lambd #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( lambd == 'loglikeratio' ) : #NEWLINE# #INDENT# a = 0 #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( lambd == 'freeman_tukey' ) : #NEWLINE# #INDENT# a = ( - 0.5 ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( lambd == 'pearson' ) : #NEWLINE# #INDENT# a = 1 #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( lambd == 'modified_loglikeratio' ) : #NEWLINE# #INDENT# a = ( - 1 ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( lambd == 'cressie_read' ) : #NEWLINE# #INDENT# a = ( 2 / 3.0 ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# raise ValueError ( ( ( 'lambd has to be a number or one of ' + 'loglikeratio, freeman_tukey, pearson, ' ) + 'modified_loglikeratio or cressie_read' ) ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# n = np . sum ( o , axis = axis ) #NEWLINE# nt = n #NEWLINE# if ( n . size > 1 ) : #NEWLINE# #INDENT# n = np . atleast_2d ( n ) #NEWLINE# if ( axis == 1 ) : #NEWLINE# #INDENT# nt = n . T #NEWLINE# #UNINDENT# if ( e . ndim == 1 ) : #NEWLINE# #INDENT# e = np . atleast_2d ( e ) #NEWLINE# if ( axis == 0 ) : #NEWLINE# #INDENT# e = e . T #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# if np . all ( ( np . sum ( e , axis = axis ) == n ) ) : #NEWLINE# #INDENT# p = ( e / ( 1.0 * nt ) ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if np . all ( ( np . sum ( e , axis = axis ) == 1 ) ) : #NEWLINE# #INDENT# p = e #NEWLINE# e = ( nt * e ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# raise ValueError ( ( 'observed and expected need to have the same ' + 'number of observations, or e needs to add to 1' ) ) #NEWLINE# #UNINDENT# #UNINDENT# k = o . shape [ axis ] #NEWLINE# if ( e . shape [ axis ] != k ) : #NEWLINE# #INDENT# raise ValueError ( ( 'observed and expected need to have the same ' + 'number of bins' ) ) #NEWLINE# #UNINDENT# if ( a == 0 ) : #NEWLINE# #INDENT# D_obs = ( ( 2 * n ) * np . sum ( ( ( o / ( 1.0 * nt ) ) * np . log ( ( o / e ) ) ) , axis = axis ) ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( a == ( - 1 ) ) : #NEWLINE# #INDENT# D_obs = ( ( 2 * n ) * np . sum ( ( ( e / ( 1.0 * nt ) ) * np . log ( ( e / o ) ) ) , axis = axis ) ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# D_obs = ( ( ( ( 2 * n ) / a ) / ( a + 1 ) ) * np . sum ( ( ( o / ( 1.0 * nt ) ) * ( ( ( o / e ) ** a ) - 1 ) ) , axis = axis ) ) #NEWLINE# #UNINDENT# #UNINDENT# return ( D_obs , stats . chi2 . sf ( D_obs , ( ( k - 1 ) - ddof ) ) )