from __future__ import division
from numpy import *
from scipy import stats
from matplotlib.pyplot import *

N = random.randn

m = 2000
n = 20

T = zeros((m,))
R = zeros((m,))
pT = zeros((m,))
pR = zeros((m,))

for k in xrange(m):
    x = N(n)
    y = N(n)

    T[k], pT[k] = stats.ttest_ind(x,y)
    R[k], pR[k] = stats.ranksums(x,y)

a = stats.t.ppf([0.025,1.-0.025], n-1)
b = stats.norm.ppf([0.025,1.-0.025])


fig = figure(figsize=(8,8),dpi=100)
ax = fig.add_axes([.3,.3,.6,.6])
axb = fig.add_axes([.3,.1,.6,.2])
axl = fig.add_axes([.1,.3,.2,.6])

ax.plot(T,R,'ok',mfc=(.7,.7,.7))
axb.hist(T,bins=50,facecolor=(1.,.7,.7),normed=True)
axl.hist(R,bins=50,facecolor=(.7,.7,1.),normed=True,orientation='horizontal')
axl.axis([0,1,-5,5])
axb.plot([a[0],a[0]],[0,1],'k--',lw=2)
axb.plot([a[1],a[1]],[0,1],'k--',lw=2)

axl.plot([0,1],[b[0],b[0]],'k--',lw=2)
axl.plot([0,1],[b[1],b[1]],'k--',lw=2)
axl.set_ylabel('standardized U statistic', fontsize=16)
axb.set_xlabel('t statistic', fontsize=16)

# print sum(1.*(T < a[0] ))/m + sum(1.*(T > a[1]))/m
# print sum(1.*(R < b[0] ))/m + sum(1.*(R > b[1]))/m

ax.fill([-5,a[0],a[0],-5],[-5,-5,5,5],color=(1.,.7,.7),alpha=.5)
ax.fill([a[1],5,5,a[1]],[-5,-5,5,5],color=(1.,.7,.7),alpha=.5)
axb.fill([-5,a[0],a[0],-5],[0,0,1,1],color=(1.,.7,.7),alpha=.5)
axb.fill([a[1],5,5,a[1]],[0,0,1,1],color=(1.,.7,.7),alpha=.5)


ax.fill([-5,-5,5,5],[-5,b[0],b[0],-5],color=(.7,.7,1.),alpha=.5)
ax.fill([-5,-5,5,5],[b[1],5,5,b[1]],color=(.7,.7,1.),alpha=.5)
axl.fill([0,0,1,1],[-5,b[0],b[0],-5],color=(.7,.7,1.),alpha=.5)
axl.fill([0,0,1,1],[b[1],5,5,b[1]],color=(.7,.7,1.),alpha=.5)



axb.axis([-5,5,0,1])
ax.axis([-5,5,-5,5])

axl.set_xticks([])
axb.set_yticks([])
axl = axl.twiny()
axb = axb.twinx()
axl.set_xticks([0,.5,1.])
axb.set_yticks([0,.5,1.])



fig.savefig('multipletesting.pdf')