Normality testing

1 files changed, 40 insertions, 0 deletions

diff --git a/data/normality.py b/data/normality.py
new file mode 100755
index 0000000..eb22553
--- /dev/null
+++ b/data/normality.py
@@ -0,0 +1,40 @@
+import sys
+import numpy as np
+from scipy.stats import lognorm
+import math
+
+def normtest(x):
+    x = np.matrix(x)
+    s = np.matrix(np.cov(x,bias=1,rowvar=0))
+    n,p = x.shape
+    dift = x - np.mean(x,axis=0)
+    dj = np.diag(dift*s.I*dift.T)
+    y = x*s.I*x.T
+    djk = - 2*y.T + np.diag(y.T).T*np.ones((1,n)) + np.ones((n,1))*np.diag(y.T)
+    b = 1/(math.sqrt(2))*((2*p + 1)/4)**(1/(p + 4))*(n**(1/(p + 4)))
+    
+    if (np.rank(s) == p):
+        hz = n * (1/(n**2) * np.sum(np.sum(np.exp( - (b**2)/2 * djk))) - 2 *
+        ((1 + (b**2))**( - p/2)) * (1/n) * (np.sum(np.exp( - ((b**2)/(2 *
+        (1 + (b**2)))) * dj))) + ((1 + (2 * (b**2)))**( - p/2)))
+    else:
+        hz = n*4
+
+    wb = (1 + b**2)*(1 + 3*b**2)
+    a = 1 + 2*b**2
+    mu = 1 - a**(- p/2)*(1 + p*b**2/a + (p*(p + 2)*(b**4))/(2*a**2))
+    si2 = 2*(1 + 4*b**2)**(- p/2) + 2*a**( - p)*(1 + (2*p*b**4)/a**2 + (3*p*
+    (p + 2)*b**8)/(4*a**4)) - 4*wb**( - p/2)*(1 + (3*p*b**4)/(2*wb) + (p*
+    (p + 2)*b**8)/(2*wb**2))
+
+    pmu = math.log(math.sqrt(mu**4/(si2 + mu**2)))
+    psi = math.sqrt(math.log((si2 + mu**2)/mu**2))
+    p = 1 - lognorm.cdf(hz,psi,scale=math.exp(pmu))
+    return p
+
+if __name__ == "__main__":
+    filename = sys.argv[1]
+    x = np.loadtxt(filename,delimiter=",")
+    x = x[:10,7:]
+    p = normtest(x)
+    print p