#BINOMIAL DISTRIBUTION
N=16
p=.7
support=0:N
pmf = dbinom(1:N,N,p)
barplot(pmf,space=0)
axis(1,0:N+.5, support)
plot(pmf)

help(hist)
n=10000
x = rbinom(n,N,p)
hist(x)
pbinom(5,N,p)
sum(dbinom(0:5,N,p))


qbinom(.7,N,p)

c(7, 9, 7, 6, 8, 6)


#BINOMIAL SAMPLES
n=50
x = rbinom(n, N, p)
plot(x)
hist(x,probability=TRUE)
boxplot(x, horizontal=TRUE)
mean(x)
Ex = sum(x)/n
median(x)

sort(x)

quantile(x)

mode(x)
var(x)
sum((x-Ex)^2)/n
sum((x-Ex)^2)/(n-1)


#NORMAL DISTRIBUTION
support = seq(40,90,.1)
mu = 65
sigma = 4
pdf = dnorm(support, mu, sigma)
plot(support,pdf)

n=16
x = rnorm(n,mu,sigma)
hist(x)

c(1, 2, 3)
ft = c(6, 5, 5, 5, 5, 5, 5,   5, 5, 5, 5, 6, 5, 5, 6, 5)
inch = c(0, 4, 6, 3, 4, 10, 10, 7, 8, 4, 8, 0, 6, 8, 2, 0)
height = 12*ft+inch
hist(height)

pnorm(60, mu, sigma)



qnorm(.5, mu, sigma)


#NORMAL SAMPLES
n=500
x = rnorm(n, mu, sigma)
x = c(x)
mean(x)
median(x)
hist(x)

plot(x)
hist(x)
boxplot(x, horizontal=TRUE)
mean(x)
Ex = sum(x)/n
median(x)
quantile(x,c(.44, .51))


mode(x)

n=50
x = rnorm(n, mu, sigma)
sqrt(var(x))
xbar = mean(x)
xbar
(sum((x-xbar)^2)/(n-1))
var(x)


sum((x-Ex)^2)/(n-1)

mu=80
sigma=8
1 - pnorm(90, mu,sigma)
pnorm(75, mu,sigma)




#TWO DATA SETS
n=50
mydata1 = rnorm(n,0,3)
mydata2 = rpois(n,4)

#THE DATA
dev.set(2)
plot(mydata2,rep(0,50))
plot(jitter(mydata2),rep(0,50))

quartz()
dev.set(3)
plot(mydata1,rep(0,50))



#MEAN/MEADIAN/ROBUSTNESS
median(mydata1)
median(mydata2)
mean(mydata1)
mean(mydata2)


bignum=NULL#-80
hist(mydata1)
lines(rep(median(c(mydata1,bignum)),2),c(0,12),lwd=3, lty="dashed")
lines(rep(mean(c(mydata1,bignum)),2),c(0,12),lwd=3)
legend("topright", c("median","mean"),lwd=c(3,3),lty=c("dashed","solid"))

bignum=NULL#40
hist(mydata2)
lines(rep(median(c(mydata2,bignum)),2),c(0,12),lwd=3, lty="dashed")
lines(rep(mean(c(mydata2,bignum)),2),c(0,12),lwd=3)
legend("topright", c("median","mean"),lwd=c(3,3),lty=c("dashed","solid"))






#For Fun...
rolls = matrix(0,6,1500)
rolls[sample(1:6,1),1] = 1

for(i in 2:1500)
{
   rolls[,i] = rolls[,i-1]
   tmp = sample(1:6,1)
   rolls[tmp,i] = rolls[tmp,i]+1
   
   plot(rolls[1,1:i]/(1:i),type='l',ylim=c(0,.5))   
   for(j in 2:6)
      lines(rolls[j,1:i]/(1:i),lty='solid')
   abline(1/6,0)
   readline()
}