Toy linefitting: bootstrapped estimator

toybootstrap.py

import numpy
from numpy import log, log10, sin, cos, tan, arctan, arccos, arcsin, abs, any, pi
import sys
import matplotlib.pyplot as plt

data = numpy.loadtxt(sys.argv[1], 
	dtype=[(colname, 'f') for colname in 'x', 'x_err', 'y', 'y_err', 'cor'],
	skiprows=1)

plt.figure(figsize=(7,7))
plt.plot(data['x'], data['y'], 'x', ms=3, color='k', alpha=0.3)
plt.xlim(data['x'].min() - 0.1, data['x'].max() + 0.1)
plt.ylim(data['y'].min() - 0.1, data['y'].max() + 0.1)
plt.xlabel('x')
plt.ylabel('y')

parameters = []
for i in range(400):
	choice = numpy.random.randint(0, len(data), size=len(data))
	sample = data[choice]
	x = numpy.mean(sample['x'])
	y = numpy.mean(sample['y'])
	x_err = numpy.std(sample['x'])
	y_err = numpy.std(sample['y'])
	corr = numpy.corrcoef([sample['x'], sample['y']])
	coef = corr[1][0]
	angle = arcsin(coef)
	# measure scatter
	deltay = (sample['y'] - y)
	deltax = (sample['x'] - x)
	xline = (deltay / sin(angle)) * cos(angle)
	distance = abs(deltax - xline) * tan(angle)
	scatter = numpy.std(distance)
	parameters.append([x, y, angle, scatter])

x_origin, y_origin, angle, scatter = numpy.median(parameters, axis=0)
logscatter = log10(scatter)
kwargs = dict(color='red')

t = numpy.linspace(-10, 10, 2)
x = t * cos(angle) + x_origin
y = t * sin(angle) + y_origin
plt.plot(x, y, '-', **kwargs)
angle2 = arctan(-1./tan(angle))
x1 = x + scatter * sin(angle2)
y1 = y + scatter * cos(angle2)
plt.plot(x1, y1, '--', **kwargs)
x1 = x - scatter * sin(angle2)
y1 = y - scatter * cos(angle2)
plt.plot(x1, y1, '--', **kwargs)
plt.plot(x_origin, y_origin, 'o', **kwargs)

plt.text(0.6, 0.8, 
	"origin: (%.2f, %.2f)\nangle: %.1f degree\nscatter: %.1f (%.1f in log)\n" % (
	x_origin, y_origin, angle * 180 / pi, scatter, logscatter),
	transform=plt.gca().transAxes)

plt.savefig('bootstrap_predict.pdf', bbox_inches='tight')
plt.savefig('bootstrap_predict.png', bbox_inches='tight')
plt.close()