from __future__ import division
import numpy as np
from scipy import stats
n = 20
pstar = 0.2
mydist = stats.binom(n,pstar)
mydist.sf(10)
k = np.arange(n+1)
mydist.sf(k)
bar(k,mydist.pmf(k),align='center');
xlim(-0.5,n+0.5)
xlabel(r'$y$');
ylabel(r'$b(y;%d,%.1f)$' % (n,pstar));
title('Binomial pmf');
savefig('notes03_binpmf.eps',bbox_inches='tight')
mydist.sf(8-1)
mypmf = mydist.pmf(k)
unlikely = mypmf < mypmf(8)
print(unlikely)
k[unlikely]
mypmf[unlikely].sum()

def ClopperPearsonCI(CL,n,x):
    tailprob = 0.5*(1.-CL)
    lower = stats.beta.ppf(tailprob,x,n-x+1)
    upper = stats.beta.isf(tailprob,x+1,n-x)
    return (lower,upper)

CIlo, CIhi = ClopperPearsonCI(.95,20,8)
CIlo
CIhi
print(stats.binom(n,CIhi).cdf(8.5))
print(stats.binom(n,CIlo).sf(7.5))

from __future__ import division
import numpy as np
from scipy import stats
n = 40
np.random.seed(1)
xi = stats.invgamma(2).rvs(size=n)
xi
p = 0.6
xpstar = 1.
teststat = np.sum(xi <= xpstar)
teststat
mu = n*p; mu
sigma = np.sqrt(n*p*(1-p)); sigma
z = (teststat-mu)/sigma; z
pvalz = 2*stats.norm.sf(z); pvalz
nulldist = stats.binom(n,p)
k = np.arange(n+1)
nullpmf = nulldist.pmf(k)
unlikely = (nullpmf <= nullpmf[teststat])
k[unlikely]
pval = np.sum(nullpmf[unlikely]); pval

n = 20
np.random.seed(1)
xi = stats.binom(4,0.4).rvs(size=n)
xi
p = 0.5
xpstar = 2
teststat1 = np.sum(xi <= xpstar); teststat1
teststat2 = np.sum(xi < xpstar); teststat2
stats.binom(n,p).sf(teststat2-0.5)
stats.binom(n,p).cdf(teststat1+0.5)
nulldist = stats.binom(n,p)
k = np.arange(n+1)
nullpmf = nulldist.pmf(k)
unlikely = (nullpmf <= nullpmf[teststat2])
np.sum(nullpmf[unlikely])
2.*stats.binom(n,p).sf(teststat2-0.5)

from __future__ import division
import numpy as np
from scipy import stats
n = 40
np.random.seed(1)
xi = stats.invgamma(2).rvs(size=n)
xi
orderstats = np.sort(xi); orderstats
p = 0.6
alpha = 1. - 0.90
mydist = stats.binom(n,p)
r,sm1 = mydist.interval(1.-alpha)
r = int(r); r
s = int(sm1) + 1; s
mu = mydist.mean(); mu
sigma = mydist.std(); sigma
zcrit = stats.norm.isf(0.5*alpha); zcrit
rn = 0.5 + mu - zcrit * sigma; rn
sn = 0.5 + mu + zcrit * sigma; sn
(int(np.floor(rn)),int(np.ceil(sn)))
(r,s)
mydist.cdf(s-1) - mydist.cdf(r-1)
mydist.cdf(s-2) - mydist.cdf(r-1)
mydist.cdf(s-1) - mydist.cdf(r)
orderstats[r-1]
orderstats[s-1]

from __future__ import division
import numpy as np
from scipy import stats
xi, yi = np.loadtxt('notes03_signtest.dat',unpack=True)
xi
yi
n = len(xi); n
jitx = stats.norm(scale=0.05).rvs(n)
jity = stats.norm(scale=0.05).rvs(n)
plot(xi+jitx,yi+jity,'k.')
xlim(0.5,6.5)
ylim(0.5,6.5)
plot([0.5,6.5],[0.5,6.5])
nplus = np.sum(xi<yi); nplus
nminus = np.sum(xi>yi); nminus
npm = nplus + nminus; npm
stats.binom(npm,0.5).cdf(nplus)
stats.binom(npm,0.5).sf(npm-nplus-0.5)
mu = 0.5*npm; mu
sigma = 0.5*np.sqrt(npm); sigma
z = (nplus - mu)/sigma; z
stats.norm.cdf(z)