p-values

Links: notebook, html, PDF, python, slides, slides(2), GitHub

Compute p-values.

from jyquickhelper import add_notebook_menu
add_notebook_menu()
%matplotlib inline

p-value table

from scipy.stats import norm
import pandas
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)
norm.ppf(0.025) -1.95996398454
-0.2 -0.1 -0.02 -0.01 -0.002 -0.001 0.001 0.002 0.01 0.02 0.1 0.2
0.001 NaN NaN NaN NaN NaN 7676 7676 1919 77 20 1.0 1.0
0.002 NaN NaN NaN NaN 3834.0 15336 15336 3834 154 39 2.0 1.0
0.050 NaN NaN 913.0 3650.0 91235.0 364939 364939 91235 3650 913 37.0 10.0
0.100 NaN 70.0 1729.0 6915.0 172866.0 691463 691463 172866 6915 1729 70.0 18.0
0.150 NaN 98.0 2449.0 9796.0 244893.0 979572 979572 244893 9796 2449 98.0 25.0
0.200 31.0 123.0 3074.0 12293.0 307317.0 1229267 1229267 307317 12293 3074 123.0 31.0
0.250 37.0 145.0 3602.0 14406.0 360137.0 1440548 1440548 360137 14406 3602 145.0 37.0
0.300 41.0 162.0 4034.0 16135.0 403354.0 1613413 1613413 403354 16135 4034 162.0 41.0
0.350 44.0 175.0 4370.0 17479.0 436966.0 1747864 1747864 436966 17479 4370 175.0 44.0
0.400 47.0 185.0 4610.0 18440.0 460976.0 1843901 1843901 460976 18440 4610 185.0 47.0
0.450 48.0 191.0 4754.0 19016.0 475381.0 1901523 1901523 475381 19016 4754 191.0 48.0
0.500 49.0 193.0 4802.0 19208.0 480183.0 1920730 1920730 480183 19208 4802 193.0 49.0
0.550 48.0 191.0 4754.0 19016.0 475381.0 1901523 1901523 475381 19016 4754 191.0 48.0
0.600 47.0 185.0 4610.0 18440.0 460976.0 1843901 1843901 460976 18440 4610 185.0 47.0
0.650 44.0 175.0 4370.0 17479.0 436966.0 1747864 1747864 436966 17479 4370 175.0 44.0
0.700 41.0 162.0 4034.0 16135.0 403354.0 1613413 1613413 403354 16135 4034 162.0 41.0
0.750 37.0 145.0 3602.0 14406.0 360137.0 1440548 1440548 360137 14406 3602 145.0 37.0
0.800 31.0 123.0 3074.0 12293.0 307317.0 1229267 1229267 307317 12293 3074 123.0 31.0
0.850 25.0 98.0 2449.0 9796.0 244893.0 979572 979572 244893 9796 2449 98.0 NaN
0.900 18.0 70.0 1729.0 6915.0 172866.0 691463 691463 172866 6915 1729 70.0 NaN
0.950 10.0 37.0 913.0 3650.0 91235.0 364939 364939 91235 3650 913 NaN NaN

p-values in 2D

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)
../_images/pvalues_examples_6_0.png ../_images/pvalues_examples_6_1.png

p-value ratio

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)
0.9513

p-values and EM

See Applying the EM Algorithm: Binomial Mixtures.

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
------- sample
ave 0.38
mea 0.364
min lk -3393.2292120130046 0.364
average pi p q likelihood
0 0.364 0.790164 0.457957 0.010195 -4146.468703
1 0.364 0.922471 0.379202 0.183127 -4040.543247
2 0.364 0.270552 0.044515 0.482497 -3468.698723
3 0.364 0.767201 0.272048 0.667035 -3245.185613
4 0.364 0.723695 0.417803 0.223080 -2977.976877
5 0.364 0.321526 0.217320 0.433511 -2899.712095
6 0.364 0.361369 0.221254 0.444773 -2843.521134
7 0.364 0.670259 0.344119 0.404412 -2825.208295
8 0.364 0.472245 0.195668 0.514627 -2823.205685
9 0.364 0.437070 0.261067 0.443920 -2752.367144

p-value and heavy tail

from scipy.stats import norm, zipf
from pyensae.mlhelper import TableFormula
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]
[6, 1, 9437, 1046824677, 1135742]
def draw_variancen(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 = TableFormula(add, columns=["index", "avg(n)", "std(n)"])
    return table.graph_XY ( [
            [ lambda v: v["index"], lambda v: v["avg(n)"], "avg(n)"],
            [ lambda v: v["index"], lambda v: v["std(n)"], "std(n)"]],
                     marker=False, link_point=True)

draw_variancen(sample)
[9999, 54761151.7614, 233970252.48746485]
<matplotlib.axes._subplots.AxesSubplot at 0x21b16175c88>
../_images/pvalues_examples_13_2.png
def draw_hill_estimator (sample) :
    res = hill_estimator(sample)
    table = TableFormula(res, columns=["d", "hill"])
    return table.graph_XY ( [
            [ lambda v: v["d"], lambda v: v["hill"], "Hill"],],
                     marker=False, link_point=True)

draw_hill_estimator(sample)
<matplotlib.axes._subplots.AxesSubplot at 0x21b12ca6470>
../_images/pvalues_examples_14_1.png
def draw_heavy_tail (sample) :
    table   = TableFormula ([ [_] for _ in sample ], columns=["obs"] )
    std = 1

    normal = norm.rvs(size = len(sample))
    normal = [ x*std  for x in normal ]
    nortbl = TableFormula ([ [_] for _ in normal ], columns=["obs"] )
    nortbl.addc ("iobs", lambda v : int(v ["obs"] * 10))
    histon  = nortbl.fgroupby ( lambda v :v["iobs"], [ lambda v : v["iobs"] ], ["nb"], [ len ] )
    histon.sort (lambda v : v["nb"], reverse = True)

    histo   = table.fgroupby ( lambda v :v["obs"], [ lambda v : v["obs"] ], ["nb"], [ len ] )

    histo.sort (lambda v : v["nb"], reverse = True)
    histo.reset_index(drop=True, inplace=True)
    histo["index"] = histo.index + 1

    vec  = list(histon["nb"])
    vec += [0,] * len(histo)
    histo.add_column_vector("nb_normal", vec [:len(histo) ] )

    histo.addc ("log(index)", lambda v : math.log(v ["index"]) / math.log(10) )
    histo.addc ("log(nb)", lambda v : math.log(v ["nb"]) / math.log(10) )
    histo.addc ("log(nb_normal)", lambda v : math.log(v ["nb_normal"]) / math.log(10) if v["nb_normal"] > 0 else 0)
    return histo.graph_XY ( [
            [ lambda v: v["log(index)"], lambda v: v["log(nb)"], "Zipf"],
            [ lambda v: v["log(index)"], lambda v: v["log(nb_normal)"], "Gauss"], ],
                     marker=False, link_point=True)

draw_heavy_tail(sample)
<matplotlib.axes._subplots.AxesSubplot at 0x21b163c7f98>
../_images/pvalues_examples_15_1.png