
import numpy as np
from scipy.stats import norm, chisqprob
from numpy.linalg import cholesky as chol
from numpy.linalg.linalg import LinAlgError
import sys
qnorm = norm.ppf
pnorm = norm.cdf

def stouffer_liptak(pvals, sigma=None, correction=False):
    """
    The stouffer_liptak correction.
    >>> stouffer_liptak([0.1, 0.2, 0.8, 0.12, 0.011])
    {'p': 0.0168..., 'C': 2.1228..., 'OK': True}

    >>> stouffer_liptak([0.5, 0.5, 0.5, 0.5, 0.5])
    {'p': 0.5, 'C': 0.0, 'OK': True}

    >>> stouffer_liptak([0.5, 0.1, 0.5, 0.5, 0.5])
    {'p': 0.28..., 'C': 0.57..., 'OK': True}

    >>> stouffer_liptak([0.5, 0.1, 0.1, 0.1, 0.5])
    {'p': 0.042..., 'C': 1.719..., 'OK': True}

    >>> stouffer_liptak([0.5], np.matrix([[1]]))
    {'p': 0.5...}
    """
    L = len(pvals)
    pvals = np.array(pvals, dtype=np.float64)
    pvals[pvals == 1] = 1.0 - 9e-16
    qvals = qnorm(1.0 - pvals, loc=0, scale=1).reshape(L, 1)
    if any(np.isinf(qvals)):
        raise Exception("bad values: %s" % pvals[list(np.isinf(qvals))])

    # dont do the correction unless sigma is specified.
    result = {"OK": True}
    if not sigma is None:
        try:
            C = chol(sigma)
            Cm1 = np.asmatrix(C).I # C^-1
            # qstar
            qvals = Cm1 * qvals
            result["OK"] = True
        except LinAlgError, e:
            print >>sys.stderr, e
            # cant do the correction non-invertible
            sigma *= 0.95
            np.fill_diagonal(sigma, 0.99)
            result["OK"] = False

        # http://en.wikipedia.org/wiki/Fisher's_method#Relation_to_Stouffer.27s_Z-score_method
    if correction:
        denom = np.sqrt(np.power(sigma, 2).sum())
        Cp = qvals.sum() / denom
    if not correction:
        Cp = qvals.sum() / np.sqrt(len(qvals))

    # get the right tail.
    pstar = norm.sf(Cp)
    if np.isnan(pstar):
        print >>sys.stderr, "BAD:", pvals, sigma
        pstar = np.median(pvals)
        result["OK"] = True
    result.update({"C": Cp, "p": pstar})
    return result

def z_score_combine(pvals, sigma):
    L = len(pvals)
    pvals = np.array(pvals, dtype=np.float64)
    pvals[pvals == 1] = 1.0 - 9e-16
    z = np.mean(qnorm(1.0 - pvals, loc=0, scale=1))
    sz = 1.0/L * np.sqrt(L + 2 * np.tril(sigma, k=-1).sum())
    res = {'p': norm.sf(z/sz), 'OK': True}
    return res

from scipy.misc import comb
from math import log, factorial
def zaykin_truncated_independent(pvals, cutoff=0.05):
    tau = cutoff
    pvals = np.asarray(pvals)
    L = len(pvals)
    # what is K?
    k = (pvals <= tau).sum()

    w = np.product(pvals[pvals <= tau]) # ? no qvalue, just a cutoff
    print w

    p = 0
    for k in range(1, L):
        lhs = comb(L, k) * (1 - tau)**(L - k)

        rhs = 0
        for s in range(k):
            rhs += ((k * log(tau) - log(w))**s) / factorial(s) * int(w <= tau**k) + tau**k * int(w > tau**k)


        p += lhs * w * rhs

    return p

def fisherp(pvals):
    """ combined fisher probability without correction """
    s = -2 * np.sum(np.log(pvals))
    return chisqprob(s, 2 * len(pvals))

if __name__ == "__main__":
    import doctest
    doctest.testmod(optionflags=doctest.ELLIPSIS)

@@ -226,26 +226,8 @@
 None
-, correction=False
 ):%0A 
@@ -834,24 +834,25 @@
 - 9e-16%0A    
+#
 qvals = qnor
@@ -896,16 +896,74 @@
 e(L, 1)%0A
+    qvals = norm.isf(pvals, loc=0, scale=1).reshape(L, 1)%0A
     if a
@@ -1304,40 +1304,8 @@
 als%0A
-            result%5B%22OK%22%5D = True%0A
     
@@ -1347,82 +1347,80 @@
     
-print %3E%3Esys.stderr, e%0A            # cant do the correction non-invertible%0A
+result%5B%22OK%22%5D = False%0A            if False:%0A                try:%0A        
     
@@ -1437,16 +1437,19 @@
 gma 
-*
+-
 = 0.
-95%0A
+03%0A   
     
@@ -1456,16 +1456,21 @@
         
+     
 np.fill_
@@ -1489,16 +1489,17 @@
 ma, 0.99
+9
 )%0A      
@@ -1508,254 +1508,426 @@
     
-result%5B%22OK%22%5D = False%0A%0A        # http://en.wikipedia.org/wiki/Fisher's_method#Relation_to_Stouffer.27s_Z-score_method%0A    if correction:%0A        denom = np.sqrt(np.power(sigma, 2).sum())%0A        Cp = qvals.sum() / denom%0A    if not correction:%0A    
+        sigma%5Bsigma %3C=0%5D = 0.001%0A                    C = chol(sigma)%0A                    Cm1 = np.asmatrix(C).I # C%5E-1%0A                    qvals = Cm1 * qvals%0A                    result%5B'OK'%5D = True%0A                except LinAlgError, e:%0A                    print %3E%3Esys.stderr, e%0A                    # cant do the correction non-invertible%0A            result = z_score_combine(pvals, sigma)%0A            return result%0A%0A
     
@@ -1965,17 +1965,16 @@
 qvals))%0A
-%0A
     # ge

use isf and try to re-test on invertible matrix