#
# Bootstrap confidence intervals
#
# Copyright (C) Wesley Phoa, 2001
#

from math import *
from Numeric import *
from RandomArray import *

def bootstrap(data, statistic, N=5000, n=None):
	"Compute statistic on N samples of size n from data."

	# if data is heavy-tailed, prefer size<len(data)
	if not n:
		n = max(12, len(data)/2)

	indices = list(map(list, randint(0, len(data), [N, n])))
	samples = map(lambda a, d=data: map(lambda i, d=d: d[i], a), indices)
	return map(statistic, samples)

def bootstrapCI(data, statistic, clevel=0.841345, N=5000, n=None):
	"Compute bootstrap confidence interval for real-valued statistic."

	results = bootstrap(data, statistic, N, n)
	try:
		results.sort()
	except:
		results = map(lambda x: x[0], results)
		results.sort()
	return(results[int(N * (1 - clevel))], results[int(N * clevel)])

def mean(data):
	"Compute mean of series."

	return sum(data) / float(len(data))

def stdev(data):
	"Compute standard deviation of series."

	l = len(data)
	d = array(data) - array(l * [mean(data)])
	return sqrt(sum(d * d) / float(l))

def correl(data):
	"Compute correlation of columns of array."

	l = len(data)
	data0 = transpose(data)[0]
	data0 = array(data0) - array(l * [mean(data0)])
	data1 = transpose(data)[1]
	data1 = array(data1) - array(l * [mean(data1)])
	var0 = sum(data0 * data0 / float(l))
	var1 = sum(data1 * data1 / float(l))
	return (sum(data0 * data1) / float(l)) / sqrt(var0 * var1)

try:
	from clipboard import theclipboard
except:
	pass

def calc(statistic):

	return statistic(theclipboard.array)

def CI(statistic, clevel=0.841345, N=5000, n=None):

	return bootstrapCI(theclipboard.array, statistic, clevel, N, n)

if __name__=='__main__':

	data = normal(0, 1, [25])
	print mean(data), bootstrapCI(data, mean)
	print stdev(data), bootstrapCI(data, stdev)

	data = normal(0, 1, [100])
	print mean(data), bootstrapCI(data, mean)
	print stdev(data), bootstrapCI(data, stdev)

	data = transpose([data, normal(0, 1, [100])])
	print correl(data), bootstrapCI(data, correl)