eval <- seq(-5,60,.001)
plot(eval,dchisq(eval,10,ncp=12), type="l", main="Population",
     xlab="Values",ylab="Density", ylim=c(0,0.06), xlim=c(0,65))


set.seed(19104)
x <- rchisq(50,10,ncp=12)
head(x)

plot(density(x), main="Sample", xlab="Value", ylim=c(0,0.06), xlim=c(0,65))

plot(density(x), main="Sample", xlab="Value", ylim=c(0,0.06), xlim=c(0,65))
points(eval,dchisq(eval,10,ncp=12),col="firebrick", type="l")

par(mfrow=c(3,3))
for(i in 1:9){
samp <- rchisq(50,10,ncp=12)
plot(density(samp), main="Sample", xlab="Value", ylim=c(0,0.06), xlim=c(0,65))
points(eval,dchisq(eval,10,ncp=12),col="firebrick", type="l")
}

par(mfrow=c(3,3))
for(i in 1:9){
samp <- rchisq(1000,10,ncp=12)
plot(density(samp), main="Sample", xlab="Value", ylim=c(0,0.06), xlim=c(0,65))
points(eval,dchisq(eval,10,ncp=12),col="firebrick", type="l")
}

mean.eval <- seq(15,29,.001)
plot(mean.eval, dnorm(mean.eval, mean=22, sd=sqrt(68/50)), type="l",
     main="Sampling Distribution of the Sample Mean", 
     xlab="Value of Sample Mean",
     ylab="Density")

sample.means <- rep(NA, 10000)

for(i in 1:10000){
  sample.means[i] <- mean(rchisq(50,10,ncp=12))
}

plot(mean.eval, dnorm(mean.eval, mean=22, sd=sqrt(68/50)), type="l",
     main="Sampling Distribution of the Sample Mean", 
     xlab="Value of Sample Mean",
     ylab="Density")
points(density(sample.means), type="l", col="darkblue")


source("https://raw.githubusercontent.com/marctrussler/IIS-Data/main/NormalShader.R")

plot(mean.eval, dnorm(mean.eval, mean=22, sd=sqrt(68/50)), type="l",
     main="Sampling Distribution of the Sample Mean", 
     xlab="Value of Sample Mean",
     ylab="Density")
normal.shader(16,28,22,sqrt(68/50),value.one = 20, value.two = 24, between=T)

pnorm(24,22,sqrt(68/50)) - pnorm(20,22,sqrt(68/50))

mean(sample.means>=20 & sample.means<=24)

qnorm(.15, 22, sqrt(68/50))

#How much above and below the mean
val <- 22 - qnorm(.15, 22, sqrt(68/50))

22 + val
22 - val



mean(sample.means>=20.79 & sample.means<=23.21)

sample.vars <- rep(NA, 10000)

for(i in 1:10000){
  sample.vars[i] <- var(rchisq(50,10,ncp=12))
}
summary(sample.vars)

library(fGarch)
set.seed(19104)
x <- rsnorm(500, mean = 73000, sd = 20000, xi = 2.2)

median(x)
mean(x)
sd(x)
var(x)

var(x)

eval <- seq(65000,75000,100)
plot(eval, dnorm(eval, mean=70000, sd=sqrt(4e+08/500)), type="l", col="firebrick")
points(eval, dnorm(eval, mean=70000, sd=sqrt(var(x)/500)), type="l", col="darkblue")


plot(eval, dnorm(eval, mean=70000, sd=sqrt(var(x)/500)), type="l", col="darkblue", ylab = "Density",
     xlab="Income")
abline(v=mean(x), lty=2, lwd=2, col="firebrick")

plot(eval, dnorm(eval, mean=70000, sd=sqrt(4e+08/500)), type="l", col="firebrick")

for(i in 1:1000){
  samp.var <- var( rsnorm(500, mean = 73000, sd = 20000, xi = 2.2))
  points(eval, dnorm(eval, mean=70000, sd=sqrt(samp.var/500)), type="l", col="gray80")
}
points(eval, dnorm(eval, mean=70000, sd=sqrt(4e+08/500)), type="l", col="firebrick")


set.seed(19104)
x <- rnorm(100,64,3)
mean(x)
var(x)

#Evaluate standard normal for alpha/2
qnorm(.3/2)

x <- seq(63,65,.001)
plot(x, dnorm(x, mean=63.92, sd=.28), type="l")
normal.shader(63,65,mu  = 63.92, sd = .28,63.62, 64.21, between=T)
abline(v=63.92)


lower <- rep(NA,10)
upper <- rep(NA,10)

for(i in 1:10){
samp <- rnorm(100,64,3)
xbar <- mean(samp)
s2 <- var(samp)
se <- sqrt(s2/100)
lower[i] <- xbar - 1.04*se
upper[i] <- xbar + 1.04*se
}

plot( seq(63.2,65,.2),seq(1,10), type="n", xlim=c(63,65))
abline(v=64, lty=2, col="firebrick")
segments(lower, 1:10, upper, 1:10)



ci.contains <- rep(NA,10000)

for(i in 1:10000){
samp <- rnorm(100,64,3)
xbar <- mean(samp)
s2 <- var(samp)
se <- sqrt(s2/100)
lower <- xbar - 1.04*se
upper<- xbar + 1.04*se
ci.contains[i] <- lower<=64 & upper>=64
}
mean(ci.contains)

qnorm(.05/2)

lower <- rep(NA,10)
upper <- rep(NA,10)

for(i in 1:10){
samp <- rnorm(100,64,3)
xbar <- mean(samp)
s <- sd(samp)
se <- s/10
lower[i] <- xbar - 1.96*se
upper[i] <- xbar + 1.96*se
}

plot( seq(63.2,65,.2),seq(1,10), type="n", xlim=c(63,65))
abline(v=64, lty=2, col="firebrick")
segments(lower, 1:10, upper, 1:10)


lower <- rep(NA,10)
upper <- rep(NA,10)

for(i in 1:10){
samp <- rnorm(10,64,3)
xbar <- mean(samp)
s <- sd(samp)
se <- s/sqrt(10)
lower[i] <- xbar - 1.04*se
upper[i] <- xbar + 1.04*se
}

plot( seq(63.2,65,.2),seq(1,10), type="n", xlim=c(63,65))
abline(v=64, lty=2, col="firebrick")
segments(lower, 1:10, upper, 1:10)



ci.contains <- rep(NA,10000)

for(i in 1:10000){
samp <- rnorm(10,64,3)
xbar <- mean(samp)
s <- sd(samp)
se <- s/sqrt(10)
lower <- xbar - 1.04*se
upper<- xbar + 1.04*se
ci.contains[i] <- lower<=64 & upper>=64
}
mean(ci.contains)

eval <- seq(-3,3,.001)
plot(eval, dnorm(eval), type="l")
points(eval, dt(eval,df=9), type="l", col="firebrick")

#From the standard normal
qnorm(.15)

#Using the t distribution
qt(.15, df=9)



ci.contains <- rep(NA,10000)

for(i in 1:10000){
samp <- rnorm(10,64,3)
xbar <- mean(samp)
s <- sd(samp)
se <- s/sqrt(10)
lower <- xbar - 1.10*se
upper<- xbar + 1.10*se
ci.contains[i] <- lower<=64 & upper>=64
}
mean(ci.contains)

eval <- seq(-3,3,.001)
plot(eval, dnorm(eval), type="l")
points(eval, dt(eval,29), type="l", col="firebrick")