chk2Dists

#!/usr/bin/python
'''
Takes two files with scores and checks if the distributions are equal
Uses the Kolmogorov-Smirnoff stat
'''
def getOptions():
    import argparse

    # create the top-level parser
    description = ("Takes two files with scores and checks if the distributions are equal. Uses the Kolmogorov-Smirnoff stat")
    parser = argparse.ArgumentParser(description = description)
    parser.add_argument('file', action='store', nargs=2,
                            help='Score file')
    parser.add_argument('-L', '--log', action="store_true",
            default=False,
            dest='log',
            help='Compute the log10 of the provided values')
    parser.add_argument('-R', '--random', action="store_true",
            default=False,
            dest='random',
            help='Compute 100 random A/B comparisons')
    parser.add_argument('-S', '--subset', action="store_true",
            default=False,
            dest='subset',
            help='Use just 1/10th of the smallest sample')
    return parser.parse_args()

options = getOptions()

import numpy as np
import sys
import math

if options.log:
	a=np.loadtxt((math.log10(float(x)) for x in open(options.file[0])))
else:
	a=np.loadtxt(open(options.file[0]))

if options.log:
        b=np.loadtxt((math.log10(float(x)) for x in open(options.file[1])))
else:
        b=np.loadtxt(open(options.file[1]))

from scipy.stats import ks_2samp

print '\t'.join(['Mean A', 'Mean B', 'StdDev A', 'StdDev B', 'p-value KS'])

i=1
import random
if options.random:
	i=100
	ps = []
if options.random or options.subset:
	if len(a) > len(b):
		sublen = len(b)/10
	else:sublen = len(a)/10

for j in range(i):
	if options.random:
		c = np.append(a,b)
		random.shuffle(c)
		if not options.subset:
			d = c[:len(a)]
			e = c[len(a):]
		else:
			d = c[:sublen]
                        e = c[sublen:sublen*2]
		k,p=ks_2samp(d, e)
		ps.append(p)
		if options.log:
			print '\t'.join([str(x) for x in [pow(10,d.mean()), pow(10,e.mean()), pow(10,d.std()), pow(10,e.std()), p]])
		else:
			print '\t'.join([str(x) for x in [d.mean(), e.mean(), d.std(), e.std(), p]])
	else:
		if options.subset:
			random.shuffle(a)
			a=a[:sublen]
			random.shuffle(b)
			b=b[:sublen]
		k,p=ks_2samp(a, b)
		if options.log:
			print '\t'.join([str(x) for x in [pow(10,a.mean()), pow(10,b.mean()), pow(10,a.std()), pow(10,b.std()), p]])
		else:
			print '\t'.join([str(x) for x in [a.mean(), b.mean(), a.std(), b.std(), p]])

if options.random:
	print '\t'.join(['Mean P-value', 'StdDev P-value'])
	print '\t'.join([str(x) for x in [np.array(ps).mean(), np.array(ps).std()]])