wrote stats.py for pval analysis

1 files changed, 38 insertions, 0 deletions

diff --git a/py/stats.py b/py/stats.py
new file mode 100644
index 0000000..b79e56f
--- /dev/null
+++ b/py/stats.py
@@ -0,0 +1,38 @@
+from data2 import trials
+from scipy.stats import cramervonmises
+from bisect import bisect_left
+from random import random
+
+def distance(pit1, pit2):
+        out = ((pit2[1]-pit1[1])**2 + (pit2[0]-pit1[0])**2)**(1/2)
+        if (out == 0):
+            return 100
+        else:
+            return out
+
+def cdf(a,nnsim,tup):
+    return [bisect_left(nnsim, x)/len(nnsim) for x in a]
+
+def calcpval(trialno):
+    size = trials[trialno].size[0]
+    pits = len(trials[trialno].pits)
+
+    nnsim = []
+    for j in range(10000):
+        curpits = []
+        for k in range(pits):
+            curpits.append([random()*size,random()*size])
+        for pit in curpits:
+            l = size
+            for pit2 in curpits:
+                l = min(l,distance(pit,pit2))
+            nnsim.append(l)
+
+    nnsim.sort()
+
+    nnreal = trials[trialno].nearest_neighbor()
+    stat = cramervonmises(nnreal,cdf,(nnsim,()))
+    return stat.pvalue
+
+for n in range(len(trials)):
+    print(calcpval(n))