from __future__ import division
import numpy as np
from scipy import stats
xi = np.array([-1.82, 0.72, 1.67, 1.09, 0.64, 0.81, 1.74, -0.80, -0.13, 1.12])
n = len(xi); n
x = np.linspace(-3,3,6001)
Fhat = np.mean(xi[None,:] <= x[:,None],axis=-1)
Fstar = stats.norm.cdf(x)
figure();
plot(x,Fhat,'b-',lw=2,label=r'$\hat{F}(x,\{x_i\})$');
plot(x,Fstar,'r--',lw=2,label=r'$F^*(x)$');
ip = np.argmax(Fstar-Fhat); x[ip]
im = np.argmax(Fhat-Fstar); x[im]
annotate('',xy=(x[ip],Fstar[ip]),xycoords='data',xytext=(x[ip],Fhat[ip]),textcoords='data',arrowprops=dict(arrowstyle='<->'));
annotate('',xy=(x[im],Fstar[im]),xycoords='data',xytext=(x[im],Fhat[im]),textcoords='data',arrowprops=dict(arrowstyle='<->'));
xlabel(r'$x$');
ylabel('cdf');
legend(loc='upper left');
Tp = max(Fstar-Fhat); Tp
Tm = max(Fhat-Fstar); Tm
x[ip+1]
xisort = np.sort(xi)
Fstari = stats.norm.cdf(xisort)
plot(xisort,Fstari,'k*',ms=10);
savefig('notes06_cdfs.eps',bbox_inches='tight');
Fhatip = np.arange(n)/n
Fhatim = (1+np.arange(n))/n
xisort[np.argmax(Fstari-Fhatip)]
max(Fstari-Fhatip)
xisort[np.argmax(Fhatim-Fstari)]
max(Fhatim-Fstari)
from scipy.special import binom
n_p = int(n*(1-Tp)); n_p
Tp * np.sum([binom(n,j)*(1-Tp-j/n)**(n-j)*(Tp+j/n)**(j-1) for j in range(n_p+1)])
n_m = int(n*(1-Tm)); n_m
Tm * np.sum([binom(n,j)*(1-Tm-j/n)**(n-j)*(Tm+j/n)**(j-1) for j in range(n_m+1)])
Tpmvals = np.linspace(1e-4,1,100)
tailprobs = np.array([Tpm * np.sum([binom(n,j)*(1-Tpm-j/n)**(n-j)*(Tpm+j/n)**(j-1) for j in range(int(n*(1-Tpm))+1)]) for Tpm in Tpmvals])
figure();
plot(Tpmvals,tailprobs);
xlabel(r'$t^{\pm}$');
ylabel(r'$P(T^{\pm}>t^{\pm})$');
title(r'Kolmogorov tail probs for $n=%d$'%n);
savefig('notes06_Ktail.eps',bbox_inches='tight');
np.exp(-2*n*Tp**2)
T95 = 0.409
T95 * np.sum([binom(n,j)*(1-T95-j/n)**(n-j)*(T95+j/n)**(j-1) for j in range(int(n*(1-T95))+1)])
F95lower = np.maximum(0.,Fhat-T95)
F95upper = np.minimum(1.,Fhat+T95)
figure();
fill_between(x,F95upper,F95lower,edgecolor='b',color=[0.9,0.9,0.9],label='95\% CI');
plot(x,Fhat,'b-',lw=2,label=r'$\hat{F}(x,\{x_i\})$');
plot(x,Fstar,'r--',lw=2,label=r'$F^*(x)$');
xlabel(r'$x$');
ylabel('cdf');
legend(loc='upper left');
savefig('notes06_CI.eps',bbox_inches='tight');

from __future__ import division
import numpy as np
from scipy import stats
from scipy.special import binom
xi = np.array([4, 1, 0, 1, 4, 3, 4])
# xi = np.array([1,1,1,2,2,2,3,3,3,3])
n = len(xi)
x = np.concatenate((range(6+1),np.arange(6+1)-1e-4))
x.sort(); x
Fxhat = np.mean(xi[None,:] <= x[:,None],axis=-1)
Fx = stats.binom(5,0.4).cdf(x)
# Fx = stats.randint(0,5).cdf(x)
Fx
figure();
plot(x,Fxhat);
plot(x,Fx);
Tp = max(Fx-Fxhat); Tp
n_p = int(n*(1-Tp)); n_p
Tp * np.sum([binom(n,j)*(1-Tp-j/n)**(n-j)*(Tp+j/n)**(j-1)
             for j in range(n_p+1)])
fj = np.array([Fx[Fx<=1-Tp-j/n][-1] for j in range(n_p)]); fj
ej = np.ones(n_p+1)
for k in (1+np.arange(n_p)): ej[k] = 1 - np.sum([binom(k,j)*fj[j]**(k-j)*ej[k] for j in range(k)])
ej
np.sum([binom(n,j)*fj[j]**(n-j)*ej[j] for j in xrange(n_p)])
Tp * np.sum([binom(n,j)*(1-Tp-j/n)**(n-j)*(Tp+j/n)**(j-1)
             for j in range(n_p+1)])

from __future__ import division
import numpy as np
from scipy import stats
xi = np.loadtxt('notes06_prelim.dat')
n = len(xi); n
xbar = np.mean(xi); xbar
s = np.std(xi); s
xi.sort()
Fhatip = np.arange(n)/n
Fhatim = (1+np.arange(n))/n
stardist = stats.norm(loc=xbar,scale=s)
Fstari = stardist.cdf(xi)
Tp = max(Fstari-Fhatip); Tp
Tm = max(Fhatim-Fstari); Tm
zi = (xi - xbar)/s
Fstarzi = stats.norm.cdf(zi)
max(Fstari-Fhatip)
max(Fhatim-Fstari)
x = np.linspace(-10,10,1000)
Fhat = np.mean(xi[None,:] <= x[:,None],axis=-1)
Fstar = stardist.cdf(x)
figure();
plot(x,Fhat,'b-',lw=2,label=r'$\hat{F}(x,\{x_i\})$');
plot(x,Fstar,'r--',lw=2,label=r'$F^*(x)$');
ip = np.argmax(Fstar-Fhat)
im = np.argmax(Fhat-Fstar)
annotate('',xy=(x[ip],Fstar[ip]),xycoords='data',xytext=(x[ip],Fhat[ip]),textcoords='data',arrowprops=dict(arrowstyle='<->'));
annotate('',xy=(x[im],Fstar[im]),xycoords='data',xytext=(x[im],Fhat[im]),textcoords='data',arrowprops=dict(arrowstyle='<->'));
xlabel(r'$x$');
ylabel('cdf');
legend(loc='upper left');
grid(True);
savefig('notes06_lilliefors.eps',bbox_inches='tight');
np.random.seed(20181030)
Nmonte = 10**5
x_Ii = stats.norm.rvs(size=(Nmonte,n))
x_Ii.sort(axis=-1)
xbar_I = np.mean(x_Ii,axis=-1)
s_I = np.std(x_Ii,axis=-1,ddof=1)
Fstar_Ii = stats.norm.cdf((x_Ii-xbar_I[:,None])/s_I[:,None])
Tp_I = np.max(Fstar_Ii-Fhatip,axis=-1)
Tm_I = np.max(Fhatim-Fstar_Ii,axis=-1)
print(np.mean(np.maximum(Tp_I,Tm_I)>max(Tp,Tm)))
xlim(-5,5);
savefig('notes06_lilliefors_zoom.eps',bbox_inches='tight');
altloc = np.median(xi); altloc
altscale = (np.percentile(xi,75) - np.percentile(xi,25))/(stats.norm.ppf(.75) - stats.norm.ppf(.25)); altscale
altdist = stats.norm(loc=altloc,scale=altscale)
Falti = altdist.cdf(xi)
altTp = max(Falti-Fhatip); altTp
altTm = max(Fhatim-Falti); altTm
Falt = altdist.cdf(x)
figure();
plot(x,Fhat,'b-',lw=2,label=r'$\hat{F}(x,\{x_i\})$');
plot(x,Falt,'r--',lw=2,label=r'$F^*(x)$');
altip = np.argmax(Falt-Fhat)
altim = np.argmax(Fhat-Falt)
annotate('',xy=(x[altip],Falt[altip]),xycoords='data',xytext=(x[altip],Fhat[altip]),textcoords='data',arrowprops=dict(arrowstyle='<->'));
annotate('',xy=(x[altim],Falt[altim]),xycoords='data',xytext=(x[altim],Fhat[altim]),textcoords='data',arrowprops=dict(arrowstyle='<->'));
xlabel(r'$x$');
ylabel('cdf');
legend(loc='upper left');
grid(True)
xlim(-5,5)
savefig('notes06_altlilliefors_zoom.eps',bbox_inches='tight')
print(np.mean(np.maximum(Tp_I,Tm_I)>max(altTp,altTm)))
TCvM = 1/(12.*n) 

from __future__ import division
import numpy as np
from scipy import stats
abfahrt = np.array([2,10,12,21,22,27,34,40,51,52,57,59])
x = np.concatenate(([abfahrt[0]],abfahrt[1:]-abfahrt[:-1])); x