Assignment_week03_June Yao_500316995

normal.data <- rnorm(100, mean = 10, sd = 10)
hist(normal.data)

par(mfrow=c(1,2)) 
set.seed(42); 
normal.data2  <-rnorm(100, mean = 10, sd = 10)
hist(normal.data2)
set.seed(42); 
normal.data3  <-rnorm(100, mean = 10, sd = 10)
hist(normal.data3)

identical(normal.data2,normal.data3)
## [1] TRUE
all.equal(normal.data2,normal.data3)
## [1] TRUE

normal.data4  <-rnorm(100, mean = 8, sd = 10)
LL <- function(mu,sigma=10) {
    R = suppressWarnings(dnorm(normal.data4, mu,sigma)) 
    return(-sum(log(R))) 
}
mu=1:20
#sapply(mu, LL)
plot(mu,sapply(mu, LL),type="l",ylab="Log likelihood")

height = read.delim("height-1.txt")
par(mfrow=c(2,2))
hist(height[,2], probability = TRUE, breaks = 5,main = "Histogram of The Heights", xlab = "Heights (m)")
hist(height[,2], probability = TRUE, breaks = 10,main = "Histogram of The Heights", xlab = "Heights (m)")
hist(height[,2], probability = TRUE, breaks = 15,main = "Histogram of The Heights", xlab = "Heights (m)")
hist(height[,2], probability = TRUE, breaks = 20,main = "Histogram of The Heights", xlab = "Heights (m)")

d1 <- density(height[,2], bw = 0.01, kernel = "gaussian") 
d2 <- density(height[,2], bw = 0.01, kernel = "epanechnikov") 
d3 <- density(height[,2], bw = 0.01, kernel = "triangular") 

hist(height[,2], freq = FALSE, breaks = 20, main = "Histogram of The Heights", xlab = "Heights (m)")
lines(d1$x, d1$y, col = "red") 
lines(d2$x, d2$y, col = "green2") 
lines(d3$x, d3$y, col = "purple") 
legend("topright",c("Gaussian","Epanechnikov", "Triangular"),
      lty = c(1, 1, 1), col = c("red", "green2", "purple"))

d4 <- density(height[,2], bw = 0.05, kernel = "gaussian") 
f <- approxfun(d4)
p <- integrate(f, 1.7,max(d4$x))
p

plot(d4, main = "Density estimation")
abline(v=1.7, col = "red")
abline(v=max(d4$x), col = "red")

library(stats4) 
mle.fit <- mle(LL, start = list(mu = 5, sigma = 8))
mle.fit
## 
## Call:
## mle(minuslogl = LL, start = list(mu = 5, sigma = 8))
## 
## Coefficients:
##       mu    sigma 
## 7.125134 8.996318

StandardNorm  <- rnorm(100, mean = 0, sd = 1)

ISE <- function(bw) {
    ds <- density(x = StandardNorm, from = -3, to = 3, n = 512,bw=bw)
    result <- suppressWarnings((ds$y - dnorm(ds$x))^2) 
    return(sum(result)) 
}
bw=seq(0.01, 1, length=10)
ISER=sapply(bw, ISE)
plot(bw,ISER,type="l",ylab="ISE")

ISERmin=min(ISER)
optimalbw=bw[which(ISER==ISERmin)]
#list(bw,ISER)
print(cat('min ISE',ISERmin,'\noptimal bandwidth:',as.double(optimalbw)))

optimalds=density(x = StandardNorm, from = -3, to = 3, n = 512,bw=optimalbw)
plot(optimalds$x,dnorm(optimalds$x),col='red',main="Density Estimation vs Normal Distribution Density",xlab="",ylab="density")
lines(optimalds$x,optimalds$y,type="l",xlab="")

bw <- bw.ucv(height[,2], lower = 0.005, upper = 0.05) 
d5 <- density(height[,2], bw = bw, kernel = "gaussian") 
plot(d5, main = "Density estimation")