#NEWLINE# def test_powerlaw_stats( ) : #NEWLINE# #INDENT# 'Test the powerlaw stats function.\n\n    This unit test is also a regression test for ticket 1548.\n\n    The exact values are:\n    mean:\n        mu = a / (a + 1)\n    variance:\n        sigma**2 = a / ((a + 2) * (a + 1) ** 2)\n    skewness:\n        One formula (see http://en.wikipedia.org/wiki/Skewness) is\n            gamma_1 = (E[X**3] - 3*mu*E[X**2] + 2*mu**3) / sigma**3\n        A short calculation shows that E[X**k] is a / (a + k), so gamma_1\n        can be implemented as\n            n = a/(a+3) - 3*(a/(a+1))*a/(a+2) + 2*(a/(a+1))**3\n            d = sqrt(a/((a+2)*(a+1)**2)) ** 3\n            gamma_1 = n/d\n        Either by simplifying, or by a direct calculation of mu_3 / sigma**3,\n        one gets the more concise formula:\n            gamma_1 = -2.0 * ((a - 1) / (a + 3)) * sqrt((a + 2) / a)\n    kurtosis: (See http://en.wikipedia.org/wiki/Kurtosis)\n        The excess kurtosis is\n            gamma_2 = mu_4 / sigma**4 - 3\n        A bit of calculus and algebra (sympy helps) shows that\n            mu_4 = 3*a*(3*a**2 - a + 2) / ((a+1)**4 * (a+2) * (a+3) * (a+4))\n        so\n            gamma_2 = 3*(3*a**2 - a + 2) * (a+2) / (a*(a+3)*(a+4)) - 3\n        which can be rearranged to\n            gamma_2 = 6 * (a**3 - a**2 - 6*a + 2) / (a*(a+3)*(a+4))\n    ' #NEWLINE# cases = [ ( 1.0 , ( 0.5 , ( 1.0 / 12 ) , 0.0 , ( - 1.2 ) ) ) , ( 2.0 , ( ( 2.0 / 3 ) , ( 2.0 / 36 ) , ( - 0.5656854249492473 ) , ( - 0.6 ) ) ) ] #NEWLINE# for ( a , exact_mvsk ) in cases : #NEWLINE# #INDENT# mvsk = stats . powerlaw . stats ( a , moments = 'mvsk' ) #NEWLINE# assert_array_almost_equal ( mvsk , exact_mvsk )