# read library

library(MASS)

# Read from the R dataset package

dim(cars)

names(cars)

summary(cars)

boxplot(cars)

boxplot(cars,las=1)     # with labels parallel or perpendicular to the axis

plot(cars, xlab = "Speed (mph)", ylab = "Stopping distance (ft)", main = "Scatter Plot",
      las=1, xlim = c(0, 25))

head(cars)              # First 6 data points

hist(cars$speed, main="Histogram of Speed", xlab="Speed")
hist(cars$dist, main="Histogram of Distance", xlab="Distance")

# Read your own data into R
Yield <- read.table("G:/teaching/Stat481/table715.dat", sep = "\t", header=TRUE)
summary(Yield)


# Webpage for Data From Textbook
# http://www.biz.uiowa.edu/faculty/jledolter/AppliedStatistics/ 

# Read from internet. Exercise 6.3-1
# Four fabrics and five runs on a wear tester, weight loss (in miligram)
WearTestr <- read.table("http://www.biz.uiowa.edu/faculty/jledolter/AppliedStatistics/R/Chapter6/Section6.3Exercise6.3-1WearTester.TXT", header=TRUE, sep = "\t")
head(WearTestr)
summary(WearTestr)

## Normal Distribution

# generate a random sample of size n from Normal dist. with mean 2 and variance 1

x.sample <- rnorm(n=100, mean=2, sd =1)
hist(x.sample)

# Density estimation
hist(x.sample, prob=TRUE)
lines(density(x.sample), col=2, lty=2)


# calculate probability P(X <= 1.6) for X~ N(2, 1)

pnorm(q=1.6, mean=2, sd=1)

# calculate quantile point P(X <=q) = 0.9 for X~ N(2, 1)

qnorm(p=0.9, mean=2, sd=1)

# density function value at x=3.1 for  X~ N(2, 1)

dnorm(x=3.1, mean=2, sd=1)

# similar for t-distribution

rt(n=10, df=5)   # random sample
pt(q=1.5, df=4)  # probability P( t(4) <= 1.5 )
qt(p=0.2, df=10) # quantile point P( t(10) <= q) =0.2
dt(x=0, df=6)    # density value at x=0 for t(6)

# Other distributions: pois(lambda), unif(min,max), f(df1, df2), chisq(df)

# Uniform discrete distribution
x = 1:20
sample(x, size=15, replace=TRUE)

xbar = 1:300
n=50

for (i in 1:300) {
  xbar[i] = mean(sample(x, size=n, replace=TRUE))
}

hist(xbar, freq=FALSE, main="Histogram of Sample Mean")
lines(density(sort(xbar)), col=2, lty=2)

# Binomial Distribution B(n,p) 

n = 100
p = 0.4

n_sim=200

x.binom <- rbinom(n_sim, n, p )
xnew <-sort(x.binom)
plot(xnew)

mu = n*p
sigma = sqrt(n*p*(1-p))

hist(xnew, prob=TRUE, main="Normal Approx. to Binoimal")
lines(density(xnew), col=2, lty=2)
lines(xnew, dnorm(xnew, mean=mu, sd=sigma), col=3, lty=1)

legend("topright", c("data","normal"), col=c(2,3), lty=c(2,1))