.. _pvaluesexamplesrst:

========
p-values
========


.. only:: html

    **Links:** :download:`notebook <pvalues_examples.ipynb>`, :downloadlink:`html <pvalues_examples2html.html>`, :download:`PDF <pvalues_examples.pdf>`, :download:`python <pvalues_examples.py>`, :downloadlink:`slides <pvalues_examples.slides.html>`, :githublink:`GitHub|_doc/notebooks/metric/pvalues_examples.ipynb|*`


Compute p-values and heavy tails estimators.

.. code:: ipython3

    from jyquickhelper import add_notebook_menu
    add_notebook_menu()






.. contents::
    :local:





.. code:: ipython3

    %matplotlib inline

p-value table
-------------

.. code:: ipython3

    from scipy.stats import norm
    import pandas
    from pandas import DataFrame
    import numpy
    
    
    def pvalue(p, q, N):
        theta   = abs(p-q)
        var     = p*(1-p) 
        bn      = (2*N)**0.5 * theta / var**0.5
        ret     = (1 - norm.cdf(bn))*2
        return ret
    
    def pvalue_N(p, q, alpha):
        theta   = abs(p-q)
        var     = p*(1-p) 
        rev     = abs(norm.ppf (alpha/2))
        N       = 2 * (rev * var**0.5 / theta)** 2
        return int(N+1)
    
    def alphatable(ps, dps, alpha):
        values = []
        for p in ps :
            row=[]
            for dp in dps :
                q = p+dp
                r = pvalue_N(p,q,alpha) if 1 >= q >= 0 else numpy.nan
                row.append (r)
            values.append (row)
        return values
        
    def dataframe(ps,dps,table):
        columns = dps
        df = pandas.DataFrame(data=table, index=ps)
        df.columns = dps
        return df
            
            
    print ("norm.ppf(0.025)",norm.ppf (0.025)) # -1.9599639845400545
    ps  = [0.001, 0.002] + [ 0.05*i for i in range (1,20) ]
    dps = [ -0.2, -0.1,  -0.02, -0.01, -0.002, -0.001,
             0.2,  0.1,   0.02,  0.01,  0.002,  0.001, ]
    dps.sort()
    t = alphatable(ps, dps, 0.05)
    dataframe (ps, dps, t)


.. parsed-literal::
    norm.ppf(0.025) -1.9599639845400545






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>-0.200</th>
          <th>-0.100</th>
          <th>-0.020</th>
          <th>-0.010</th>
          <th>-0.002</th>
          <th>-0.001</th>
          <th>0.001</th>
          <th>0.002</th>
          <th>0.010</th>
          <th>0.020</th>
          <th>0.100</th>
          <th>0.200</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0.001</th>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>7676</td>
          <td>7676</td>
          <td>1919</td>
          <td>77</td>
          <td>20</td>
          <td>1.0</td>
          <td>1.0</td>
        </tr>
        <tr>
          <th>0.002</th>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>3834.0</td>
          <td>15336</td>
          <td>15336</td>
          <td>3834</td>
          <td>154</td>
          <td>39</td>
          <td>2.0</td>
          <td>1.0</td>
        </tr>
        <tr>
          <th>0.050</th>
          <td>NaN</td>
          <td>NaN</td>
          <td>913.0</td>
          <td>3650.0</td>
          <td>91235.0</td>
          <td>364939</td>
          <td>364939</td>
          <td>91235</td>
          <td>3650</td>
          <td>913</td>
          <td>37.0</td>
          <td>10.0</td>
        </tr>
        <tr>
          <th>0.100</th>
          <td>NaN</td>
          <td>70.0</td>
          <td>1729.0</td>
          <td>6915.0</td>
          <td>172866.0</td>
          <td>691463</td>
          <td>691463</td>
          <td>172866</td>
          <td>6915</td>
          <td>1729</td>
          <td>70.0</td>
          <td>18.0</td>
        </tr>
        <tr>
          <th>0.150</th>
          <td>NaN</td>
          <td>98.0</td>
          <td>2449.0</td>
          <td>9796.0</td>
          <td>244893.0</td>
          <td>979572</td>
          <td>979572</td>
          <td>244893</td>
          <td>9796</td>
          <td>2449</td>
          <td>98.0</td>
          <td>25.0</td>
        </tr>
        <tr>
          <th>0.200</th>
          <td>31.0</td>
          <td>123.0</td>
          <td>3074.0</td>
          <td>12293.0</td>
          <td>307317.0</td>
          <td>1229267</td>
          <td>1229267</td>
          <td>307317</td>
          <td>12293</td>
          <td>3074</td>
          <td>123.0</td>
          <td>31.0</td>
        </tr>
        <tr>
          <th>0.250</th>
          <td>37.0</td>
          <td>145.0</td>
          <td>3602.0</td>
          <td>14406.0</td>
          <td>360137.0</td>
          <td>1440548</td>
          <td>1440548</td>
          <td>360137</td>
          <td>14406</td>
          <td>3602</td>
          <td>145.0</td>
          <td>37.0</td>
        </tr>
        <tr>
          <th>0.300</th>
          <td>41.0</td>
          <td>162.0</td>
          <td>4034.0</td>
          <td>16135.0</td>
          <td>403354.0</td>
          <td>1613413</td>
          <td>1613413</td>
          <td>403354</td>
          <td>16135</td>
          <td>4034</td>
          <td>162.0</td>
          <td>41.0</td>
        </tr>
        <tr>
          <th>0.350</th>
          <td>44.0</td>
          <td>175.0</td>
          <td>4370.0</td>
          <td>17479.0</td>
          <td>436966.0</td>
          <td>1747864</td>
          <td>1747864</td>
          <td>436966</td>
          <td>17479</td>
          <td>4370</td>
          <td>175.0</td>
          <td>44.0</td>
        </tr>
        <tr>
          <th>0.400</th>
          <td>47.0</td>
          <td>185.0</td>
          <td>4610.0</td>
          <td>18440.0</td>
          <td>460976.0</td>
          <td>1843901</td>
          <td>1843901</td>
          <td>460976</td>
          <td>18440</td>
          <td>4610</td>
          <td>185.0</td>
          <td>47.0</td>
        </tr>
        <tr>
          <th>0.450</th>
          <td>48.0</td>
          <td>191.0</td>
          <td>4754.0</td>
          <td>19016.0</td>
          <td>475381.0</td>
          <td>1901523</td>
          <td>1901523</td>
          <td>475381</td>
          <td>19016</td>
          <td>4754</td>
          <td>191.0</td>
          <td>48.0</td>
        </tr>
        <tr>
          <th>0.500</th>
          <td>49.0</td>
          <td>193.0</td>
          <td>4802.0</td>
          <td>19208.0</td>
          <td>480183.0</td>
          <td>1920730</td>
          <td>1920730</td>
          <td>480183</td>
          <td>19208</td>
          <td>4802</td>
          <td>193.0</td>
          <td>49.0</td>
        </tr>
        <tr>
          <th>0.550</th>
          <td>48.0</td>
          <td>191.0</td>
          <td>4754.0</td>
          <td>19016.0</td>
          <td>475381.0</td>
          <td>1901523</td>
          <td>1901523</td>
          <td>475381</td>
          <td>19016</td>
          <td>4754</td>
          <td>191.0</td>
          <td>48.0</td>
        </tr>
        <tr>
          <th>0.600</th>
          <td>47.0</td>
          <td>185.0</td>
          <td>4610.0</td>
          <td>18440.0</td>
          <td>460976.0</td>
          <td>1843901</td>
          <td>1843901</td>
          <td>460976</td>
          <td>18440</td>
          <td>4610</td>
          <td>185.0</td>
          <td>47.0</td>
        </tr>
        <tr>
          <th>0.650</th>
          <td>44.0</td>
          <td>175.0</td>
          <td>4370.0</td>
          <td>17479.0</td>
          <td>436966.0</td>
          <td>1747864</td>
          <td>1747864</td>
          <td>436966</td>
          <td>17479</td>
          <td>4370</td>
          <td>175.0</td>
          <td>44.0</td>
        </tr>
        <tr>
          <th>0.700</th>
          <td>41.0</td>
          <td>162.0</td>
          <td>4034.0</td>
          <td>16135.0</td>
          <td>403354.0</td>
          <td>1613413</td>
          <td>1613413</td>
          <td>403354</td>
          <td>16135</td>
          <td>4034</td>
          <td>162.0</td>
          <td>41.0</td>
        </tr>
        <tr>
          <th>0.750</th>
          <td>37.0</td>
          <td>145.0</td>
          <td>3602.0</td>
          <td>14406.0</td>
          <td>360137.0</td>
          <td>1440548</td>
          <td>1440548</td>
          <td>360137</td>
          <td>14406</td>
          <td>3602</td>
          <td>145.0</td>
          <td>37.0</td>
        </tr>
        <tr>
          <th>0.800</th>
          <td>31.0</td>
          <td>123.0</td>
          <td>3074.0</td>
          <td>12293.0</td>
          <td>307317.0</td>
          <td>1229267</td>
          <td>1229267</td>
          <td>307317</td>
          <td>12293</td>
          <td>3074</td>
          <td>123.0</td>
          <td>31.0</td>
        </tr>
        <tr>
          <th>0.850</th>
          <td>25.0</td>
          <td>98.0</td>
          <td>2449.0</td>
          <td>9796.0</td>
          <td>244893.0</td>
          <td>979572</td>
          <td>979572</td>
          <td>244893</td>
          <td>9796</td>
          <td>2449</td>
          <td>98.0</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>0.900</th>
          <td>18.0</td>
          <td>70.0</td>
          <td>1729.0</td>
          <td>6915.0</td>
          <td>172866.0</td>
          <td>691463</td>
          <td>691463</td>
          <td>172866</td>
          <td>6915</td>
          <td>1729</td>
          <td>70.0</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>0.950</th>
          <td>10.0</td>
          <td>37.0</td>
          <td>913.0</td>
          <td>3650.0</td>
          <td>91235.0</td>
          <td>364939</td>
          <td>364939</td>
          <td>91235</td>
          <td>3650</td>
          <td>913</td>
          <td>NaN</td>
          <td>NaN</td>
        </tr>
      </tbody>
    </table>
    </div>



p-values in 2D
--------------

.. code:: ipython3

    import numpy, matplotlib, random, math
    import matplotlib.pyplot as pylab
    
    
    def matrix_square_root(sigma) :
        eigen, vect = numpy.linalg.eig(sigma)
        dim = len(sigma)
        res = numpy.identity(dim)
        for i in range(0,dim) :
            res[i,i] = eigen[i]**0.5
        return vect * res * vect.transpose()
        
    def chi2_level (alpha = 0.95) :
        N = 1000
        x = [ random.gauss(0,1) for _ in range(0,N) ]
        y = [ random.gauss(0,1) for _ in range(0,N) ]
        r = map ( lambda c : (c[0]**2+c[1]**2)**0.5, zip(x,y))
        r = list(r)
        r.sort()
        res = r [ int (alpha * N) ]
        return res
        
    def square_figure(mat, a) :    
        x = [ ]
        y = [ ]
        for i in range (0,100) :
            x.append( a * mat[0][0]**0.5 ) 
            y.append( (random.random ()-0.5) * a * mat[1][1]**0.5*2 )
            x.append( -a * mat[0][0]**0.5 ) 
            y.append( (random.random ()-0.5) * a * mat[1][1]**0.5*2 )
    
            y.append( a * mat[1][1]**0.5 ) 
            x.append( (random.random ()-0.5) * a * mat[0][0]**0.5*2 )
            y.append( -a * mat[1][1]**0.5 ) 
            x.append( (random.random ()-0.5) * a * mat[0][0]**0.5*2 )
            
        pylab.plot(x,y, 'ro')
        
        x = [ ]
        y = [ ]
        for i in range (0,100) :
            x.append( a ) 
            y.append( (random.random ()-0.5) * a*2 )
            x.append( -a ) 
            y.append( (random.random ()-0.5) * a*2  )
            
            y.append( a ) 
            x.append( (random.random ()-0.5) * a*2  )
            y.append( -a ) 
            x.append( (random.random ()-0.5) * a*2  )
            
        xs,ys = [],[]
        for a,b in zip (x,y) :
            ar = numpy.matrix( [ [a], [b] ] ).transpose()
            we = ar * root
            xs.append( we [0,0] )
            ys.append( we [0,1] )
            
        pylab.plot(xs,ys, 'bo')
        pylab.show()
    
    def circle_figure (mat, a) :
        x = [ ]
        y = [ ]
        for i in range (0,200) :
            z = random.random() * math.pi * 2
            i = a * mat[0][0]**0.5 * math.cos(z)
            j = a * mat[0][0]**0.5 * math.sin(z)
            x.append ( i )
            y.append ( j )
        pylab.plot(x,y, 'ro')
        
        x = [ ]
        y = [ ]
        for i in range (0,200) :
            z = random.random() * math.pi * 2
            i = a * math.cos(z)
            j = a * math.sin(z)
            x.append ( i )
            y.append ( j )
            
        xs,ys = [],[]
        for a,b in zip (x,y) :
            ar = numpy.matrix( [ [a], [b] ] ).transpose()
            we = ar * root
            xs.append( we [0,0] )
            ys.append( we [0,1] )
            
        pylab.plot(xs,ys, 'bo')
        pylab.show()
    
    level = chi2_level ()
    mat   = [ [0.1, 0.05], [0.05, 0.2] ]
    npmat = numpy.matrix(mat)
    root  = matrix_square_root (npmat)
    square_figure (mat, 1.96)
    circle_figure (mat, level)



.. image:: pvalues_examples_6_0.png



.. image:: pvalues_examples_6_1.png


p-value ratio
-------------

.. code:: ipython3

    import random, math
    
    def densite_gauss (mu, sigma, x) :
        e = -(x - mu)**2 / (sigma**2 * 2)
        d = 1. / ((2*math.pi)**0.5 * sigma)
        return d * math.exp (e)
        
    def simulation_vector (N, mu, sigma) :
        return [ random.gauss(mu,sigma) for n in range(N) ]
        
    def ratio (vector, x, fdensite) :
        under = 0
        above = 0
        fx    = fdensite(x)
        for u in vector:
            f = fdensite (u)
            if f >= fx:
                above += 1
            else:
                under += 1
        return float(above) / float (above + under)
            
    x     = 1.96
    N     = 10000
    mu    = 0
    sigma = 1
    
    v = simulation_vector(N, mu, sigma)
    g = ratio(v, x, lambda y: densite_gauss (mu, sigma, y) )
    print (g)


.. parsed-literal::
    0.9487


p-values and EM
---------------

See `Applying the EM Algorithm: Binomial
Mixtures <http://statisticalrecipes.blogspot.fr/2012/04/applying-em-algorithm-binomial-mixtures.html>`__.

.. code:: ipython3

    from scipy.stats import norm
    import random, math
    
    
    def average_std_deviation(sample):
        mean = 0.
        var  = 0.
        for x in sample:
            mean += x
            var += x*x
        mean /= len(sample)
        var /= len(sample)
        var -= mean*mean
        return mean,var ** 0.5
    
    def bootsample(sample):
        n = len(sample)-1
        return [ sample[ random.randint(0,n) ] for _ in sample ]
        
    def bootstrap_difference(sampleX, sampleY, draws=2000, confidence=0.05):
        diff = [ ]
        for n in range (0,draws) :
            if n % 1000 == 0: 
                print(n)
            sx = bootsample(sampleX)
            sy = bootsample(sampleY)
            px = sum(sx) * 1.0/ len(sx)
            py = sum(sy) * 1.0/ len(sy)
            diff.append (px-py) 
        diff.sort()
        n = int(len(diff) * confidence / 2)
        av = sum(diff) / len(diff)
        return av, diff [n], diff [len(diff)-n]
    
    # generation of a sample
    
    def generate_obs(p):
        x = random.random()
        if x <= p : return 1
        else : return 0
    
    def generate_n_obs(p, n):
        return [ generate_obs(p) for i in range (0,n) ]
            
    # std deviation
    
    def diff_std_deviation(px, py):
        s = px*(1-px) + py*(1-py)
        return px, py, s**0.5
    
    def pvalue(diff, std, N):
        theta = abs(diff)
        bn = (2*N)**0.5 * theta / std
        pv = (1 - norm.cdf(bn))*2
        return pv
        
    def omega_i (X, pi, p, q) :
        np = p * pi     if X == 1 else (1-p)*pi
        nq = q * (1-pi) if X == 1 else (1-q)*(1-pi)
        return np / (np + nq)
        
    def likelihood (X, pi, p, q) :
        np = p * pi     if X == 1 else (1-p)*pi
        nq = q * (1-pi) if X == 1 else (1-q)*(1-pi)
        return math.log(np) + math.log(nq)
        
    def algoEM (sample):
        p   = random.random()
        q   = random.random()
        pi  = random.random()
        iter   = 0
        while iter < 10 :
            lk = sum ( [ likelihood (x, pi, p, q) for x in sample ] )
            wi  = [ omega_i (x, pi, p, q) for x in sample ]
            sw  = sum(wi)
            pin = sum(wi) / len(wi)
            pn  = sum([ x * w     for x,w in zip (sample,wi) ]) / sw
            qn  = sum([ x * (1-w) for x,w in zip (sample,wi) ]) / (len(wi) - sw)
            
            pi,p,q = pin,pn,qn
            iter += 1
            
        lk = sum ( [ likelihood (x, pi, p, q) for x in sample ] )
        return pi,p,q, lk
    
    
    # mix
    p,q   = 0.20, 0.80
    pi    = 0.7
    N     = 1000
    na    = int(N * pi)
    nb    = N - na
    
    print("------- sample")
    sampleX = generate_n_obs(p, na) + generate_n_obs (q, nb)
    random.shuffle(sampleX)
    print("ave", p * pi + q*(1-pi))
    print("mea", sum(sampleX)*1./len(sampleX))
    
    lk = sum ( [ likelihood (x, pi, p, q) for x in sampleX ] )
    print ("min lk", lk, sum (sampleX)*1. / len(sampleX))
    res = []
    for k in range (0, 10) :
        r = algoEM (sampleX)
        res.append ( (r[-1], r) )
    res.sort ()
    
    rows = []
    for r in res:
        pi,p,q,lk = r[1]
        rows.append( [p * pi + q*(1-pi)] + list(r[1]))
    
    df = pandas.DataFrame(data=rows)
    df.columns = ["average", "pi", "p", "q", "likelihood"]
    df


.. parsed-literal::
    ------- sample
    ave 0.38
    mea 0.373
    min lk -3393.2292120130046 0.373






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>average</th>
          <th>pi</th>
          <th>p</th>
          <th>q</th>
          <th>likelihood</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>0.373</td>
          <td>0.000324</td>
          <td>0.341877</td>
          <td>0.373010</td>
          <td>-9358.705695</td>
        </tr>
        <tr>
          <th>1</th>
          <td>0.373</td>
          <td>0.863747</td>
          <td>0.284788</td>
          <td>0.932204</td>
          <td>-4531.967709</td>
        </tr>
        <tr>
          <th>2</th>
          <td>0.373</td>
          <td>0.936083</td>
          <td>0.346101</td>
          <td>0.766941</td>
          <td>-4490.512057</td>
        </tr>
        <tr>
          <th>3</th>
          <td>0.373</td>
          <td>0.123023</td>
          <td>0.290964</td>
          <td>0.384508</td>
          <td>-3563.557269</td>
        </tr>
        <tr>
          <th>4</th>
          <td>0.373</td>
          <td>0.538835</td>
          <td>0.053584</td>
          <td>0.746213</td>
          <td>-3487.438442</td>
        </tr>
        <tr>
          <th>5</th>
          <td>0.373</td>
          <td>0.346351</td>
          <td>0.057880</td>
          <td>0.539974</td>
          <td>-3302.391944</td>
        </tr>
        <tr>
          <th>6</th>
          <td>0.373</td>
          <td>0.797540</td>
          <td>0.376491</td>
          <td>0.359248</td>
          <td>-3144.938682</td>
        </tr>
        <tr>
          <th>7</th>
          <td>0.373</td>
          <td>0.392520</td>
          <td>0.592563</td>
          <td>0.231131</td>
          <td>-2902.915478</td>
        </tr>
        <tr>
          <th>8</th>
          <td>0.373</td>
          <td>0.390241</td>
          <td>0.459488</td>
          <td>0.317648</td>
          <td>-2778.903072</td>
        </tr>
        <tr>
          <th>9</th>
          <td>0.373</td>
          <td>0.609127</td>
          <td>0.338062</td>
          <td>0.427447</td>
          <td>-2764.987703</td>
        </tr>
      </tbody>
    </table>
    </div>



p-value and heavy tail
----------------------

.. code:: ipython3

    from scipy.stats import norm, zipf
    import sys
    
    
    def generate_n_obs_zipf (tail_index, n) :
        return list(zipf.rvs(tail_index, size=n))    
                         
    def hill_estimator (sample) :
        sample = list(sample)
        sample.sort(reverse=True)
        end = len(sample)/10
        end = min(end,100)
        s = 0.
        res = []
        for k in range (0,end) :
            s += math.log(sample[k])
            h = (s - (k+1)*math.log(sample[k+1]))/(k+1)
            h = 1./h
            res.append( [k, h] )
        return res
        
    # mix
    tail_index  = 1.05
    N           = 10000
    
    sample = generate_n_obs_zipf(tail_index, N)
    sample[:5]




.. parsed-literal::
    [357621, 148, 18, 1812876449, 36150]



.. code:: ipython3

    import pandas
    
    
    def graph_XY(curves, xlabel=None, ylabel=None, marker=True,
                 link_point=False, title=None, format_date="%Y-%m-%d",
                 legend_loc=0, figsize=None, ax=None):
        if ax is None:
            import matplotlib.pyplot as plt  # pylint: disable=C0415
            fig, ax = plt.subplots(1, 1, figsize=figsize)
    
        smarker = {(True, True): 'o-', (True, False): 'o', (False, True): '-',
                   # (False, False) :''
                   }[marker, link_point]
        has_date = False
        for xf, yf, label in curves:
            ax.plot(xf, yf, smarker, label=label)
        ax.legend(loc=legend_loc)
        return ax
    
    
    def draw_variance(sample) :
        avg = 0.
        std = 0.
        n   = 0.
        w   = 1.
        add = [] 
        for i,x in enumerate(sample) :
            x    = float (x)
            avg += x * w
            std += x*x * w
            n   += w
            val  = (std/n - (avg/n)**2)**0.5
            add.append ( [ i, avg/n, val] )
        
        print(add[-1])
        table = pandas.DataFrame(add, columns=["index", "avg(n)", "std(n)"])
        return graph_XY([
                [table['index'], table["avg(n)"], "avg(n)"],
                [table['index'], table["std(n)"], "std(n)"],
            ], marker=False, link_point=True)  
    
    draw_variance(sample);


.. parsed-literal::
    [9999, 55186871.0339, 233342554.46156308]



.. image:: pvalues_examples_13_1.png


.. code:: ipython3

    def draw_hill_estimator (sample) :
        res = hill_estimator(sample)
        table = DataFrame(res, columns=["d", "hill"])
        return graph_XY(
            [[table['d'], table['hill'], "Hill"],],
            marker=False, link_point=True)
    
    draw_hill_estimator(sample);



.. image:: pvalues_examples_14_0.png


.. code:: ipython3

    def draw_heavy_tail (sample) :
        table   = DataFrame([[_] for _ in sample ], columns=["obs"])
        std = 1
    
        normal = norm.rvs(size = len(sample))
        normal = [ x*std  for x in normal ]
        nortbl = DataFrame([ [_] for _ in normal ], columns=["obs"])
        nortbl["iobs"] = (nortbl['obs'] * 10).astype(numpy.int64)
    
        histon = nortbl[["iobs", "obs"]].groupby('iobs', as_index=False).count()
        histon.columns = ['iobs', 'nb']
        histon = histon.sort_values("nb", ascending=False).reset_index(drop=True)
        
        table["one"] = 1
        histo = table.groupby('obs', as_index=False).count()
        histo.columns = ['obs', 'nb']                            
        histo = histo.sort_values('nb', ascending=False).reset_index(drop=True)
        histo.reset_index(drop=True, inplace=True)
        histo["index"] = histo.index + 1
        
        vec  = list(histon["nb"])
        vec += [0,] * len(histo)
        histo['nb_normal'] = vec[:len(histo)]
        
        histo["log(index)"] = numpy.log(histo["index"]) / numpy.log(10)
        histo["log(nb)"] = numpy.log(histo["nb"]) / numpy.log(10)
        histo["log(nb_normal)"] = numpy.log(histo["nb_normal"]) / numpy.log(10)
        return graph_XY ([
                [histo["log(index)"], histo["log(nb)"], "Zipf"],
                [histo["log(index)"], histo["log(nb_normal)"], "Gauss"], ],
                         marker=False, link_point=True)  
    
    draw_heavy_tail(sample);


.. parsed-literal::
    c:\python372_x64\lib\site-packages\pandas\core\series.py:679: RuntimeWarning: divide by zero encountered in log
      result = getattr(ufunc, method)(*inputs, **kwargs)



.. image:: pvalues_examples_15_1.png