# Unit 1 for R4STAT				030312
#				prepared by C. Andy Tsao
# Using R for data exploration
data(faithful)     # Load old faithful data
summary(faithful)  # Brief summary
attach(faithful)   # Use the data set as default
stem(eruptions)
hist(eruptions)
hist(eruptions,prob=T)
lines(density(eruptions))
detach(faithful)

# Calculating mean and variance for discrete r.v.
# Ex. 5.37 of Johnson and Bhattcharyyya

pr <- c(.315,.289,.201,.114,.063,.012,.006)
x <-seq(0,6,1)
pr
x
sum(pr)       # Check pr is a prob assignment
mean(x)       # plain average without prob
ex<-x%*%pr;ex    # expectation, E(X)
v<-((x-ex)^2)%*%pr;v # Var(X)
sdv<-sqrt(v); sdv    # Standard Dev
# Alternative formula for
ex2<-(x^2)%*%pr      # Var(X)= E(X^2) - (EX)^2
ex2-ex^2
v

rm(list=ls(all=TRUE)) # Clear up the variables

# Some experiment for binomial r.v.
# How does it look like (pmf, cdf)? X \sim Bin(n,p)

n<-20;p<-0.5; q<-1-p;
x<-seq(0,n,1)
par(mfrow=c(2,2)); # Set up graphical configuration
plot(dbinom(x,n,p)~x,type="h")  # P(X = x)
plot(pbinom(x,n,p)~x,type="h")  # P(X \leq x)
windows()				    # Open a new window
plot(dbinom(x,n,p)~x,type="h")  # P(X = x)
plot(pbinom(x,n,p)~x,type="h")  # P(X \leq x)
ex<-x%*%dbinom(x,n,p);ex        # Checking EX=np
                                # Var(X) = npq
vx<-((x-ex)^2)%*%dbinom(x,n,p);vx
n*p
n*p*q
# Generate 15 random numbers from Bin(n,p)
plot(dbinom(x,n,p)~x,type="h")  # P(X = x)
r30<-rbinom(30,n,p);
r60<-rbinom(60,n,p);
r90<-rbinom(90,n,p);
hist(r30,xlim=c(0,n),prob=T);
hist(r60,xlim=c(0,n),prob=T);
hist(r90,xlim=c(0,n),prob=T);
s30<-(r30-n*p)/sqrt(n*p*q)
s60<-(r60-n*p)/sqrt(n*p*q)
s90<-(r90-n*p)/sqrt(n*p*q)
hist(s30,prob=T);
hist(s60,prob=T);
hist(s90,prob=T);
# Very large number
r10000<-rbinom(10000,n,p);
s10000<-(r10000-n*p)/sqrt(n*p*q) ;
hist(s10000,prob=T);

# For those who are into lottery
no<-seq(1,42,1)
sort(sample(no,6))

# Have fun!