GaryBAYLOR / R-code

A collection of algorithms written by myself for solving statistical problems
0 stars 0 forks source link

Coding is funny, Coding is thrilling, Coding is great ! #6

Open GaryBAYLOR opened 9 years ago

GaryBAYLOR commented 9 years ago

The more I code using R, the more I found it amazing. Sometimes just an idea comes to my mind and I wonder can I fulfil it using R? I am not sure if I can succeed and I just have an picture in my mind. After I spend an hour or a lot more time on coding the great moment come ! R gives me the same output as I want, and even more. So how can I not love R? The following is a tiny piece of code I wrote this afternoon to visualize the coverage of confidence interval of MLE. n = 20 image

n = 50 image

n=100 image

Attach the code.

coverage <- function(repetetion = 20, n = 11, mu = 10, sigma = 3) {
    lower <- upper <- include <- numeric(repetetion)
    for(i in 1:repetetion) {
        x <- rnorm(n, mu, sigma)
        meanx <- mean(x)
        sigma2x <- sum((x - meanx)^2)/n
        lower[i] <- meanx - sd(x) * qt(1 - 0.05/2, n-1) / sqrt(n)
        upper[i] <- meanx + sd(x) * qt(1 - 0.05/2, n-1) / sqrt(n)
        if(lower[i] <= mu & upper[i] >= mu) {
            include[i] <- 1
        }
        else {
            include[i] <- 0
        }       
    }
    maxx <- max(upper)
    minx <- min(lower)
    plot(0, 0, xlim = c(0, repetetion+1), ylim = c(minx, maxx), type = "n", xlab = "repetetion", ylab = "Coverage", main = "Confidence Interval Coverage")
    for(i in 1:repetetion) {
        xi <- rep(i, 10)
        yi <- seq(from = lower[i], to = upper[i], length = 10)
        lines(xi, yi, lty = 1, lwd = 1.5, col = ifelse(include[i] == 0, "red", "black"))
    }
    return(list(total = repetetion, covered = sum(include), coverage = sum(include) / repetetion))
}