138 lines
3.6 KiB
R
138 lines
3.6 KiB
R
load(url("http://ls.home.amu.edu.pl/data_sets/Centrala.RData"))
|
|
library(EnvStats)
|
|
epois(Centrala$Liczba,
|
|
method = "mle/mme/mvue",
|
|
ci = TRUE, ci.type = "two-sided", conf.level = 0.95,
|
|
ci.method = "exact")$interval$limits
|
|
epois(Centrala$Liczba,
|
|
method = "mle/mme/mvue",
|
|
ci = TRUE, ci.type = "two-sided", conf.level = 0.95,
|
|
ci.method = "pearson.hartley.approx")$interval$limits
|
|
epois(Centrala$Liczba,
|
|
method = "mle/mme/mvue",
|
|
ci = TRUE, ci.type = "two-sided", conf.level = 0.95,
|
|
ci.method = "normal.approx")$interval$limits
|
|
View(Centrala)
|
|
awarie <- read.table("http://ls.home.amu.edu.pl/data_sets/awarie.txt")
|
|
lambda_limits <- eexp(awarie$V1,
|
|
method = "mle/mme",
|
|
ci = TRUE, ci.type = "two-sided", conf.level = 0.95,
|
|
ci.method = "exact")$interval$limits
|
|
rev(1 / lambda_limits)
|
|
rev(1 / lambda_limits^2)
|
|
median_cint <- function(x, conf_level = 0.95) {
|
|
m <- mean(x^2)
|
|
a <- 1 - conf_level
|
|
n <- length(x)
|
|
L <- sqrt(log(2) * m * (1 - qnorm(1 - a / 2) / sqrt(n)))
|
|
R <- sqrt(log(2) * m * (1 + qnorm(1 - a / 2) / sqrt(n)))
|
|
w <- list(title = "mediana", est = sqrt(log(2) * m), l = L, r = R, conf_level = conf_level)
|
|
class(w) <- "confint"
|
|
return(w)
|
|
}
|
|
print.confint <- function(x) {
|
|
cat(x$conf_level * 100, "percent confidence interval:", "\n")
|
|
cat(x$l, " ", x$r, "\n")
|
|
}
|
|
summary.confint <- function(x) {
|
|
cat("\n", "Confidence interval of", x$title, "\n", "\n")
|
|
cat(x$conf_level * 100, "percent confidence interval:", "\n")
|
|
cat(x$l, " ", x$r, "\n")
|
|
cat("sample estimate", "\n")
|
|
cat(x$est, "\n")
|
|
}
|
|
x <- c(0.9, 6.2, 2.1, 4.1, 7.3,
|
|
1.0, 4.6, 6.4, 3.8, 5.0,
|
|
2.7, 9.2, 5.9, 7.4, 3.0,
|
|
4.9, 8.2, 5.0, 1.2, 10.1,
|
|
12.2, 2.8, 5.9, 8.2, 0.5)
|
|
median_cint(x)
|
|
summary(median_cint(x))
|
|
conf_int <- function(x, conf_level = 0.95) {
|
|
c(mean(x) - sd(x) * qt(1 - (1 - conf_level) / 2, length(x) - 1) / sqrt(length(x)),
|
|
mean(x) + sd(x) * qt(1 - (1 - conf_level) / 2, length(x) - 1) / sqrt(length(x)))
|
|
}
|
|
mu <- 1
|
|
sigma <- 3
|
|
nr <- 1000
|
|
df <- 3
|
|
lambda <- 3
|
|
n <- 10
|
|
cat(paste("n = ", n, sep = ""), "\n")
|
|
temp <- numeric(nr)
|
|
set.seed(1234)
|
|
for (i_nr in seq_len(nr)) {
|
|
x <- rnorm(n, mu, sigma)
|
|
temp_int <- conf_int(x)
|
|
temp[i_nr] <- ((temp_int[1] < mu) & (mu < temp_int[2]))
|
|
}
|
|
mean(temp)
|
|
temp <- numeric(nr)
|
|
set.seed(1234)
|
|
for (i_nr in seq_len(nr)) {
|
|
x <- rchisq(n, df)
|
|
temp_int <- conf_int(x)
|
|
temp[i_nr] <- ((temp_int[1] < df) & (df < temp_int[2]))
|
|
}
|
|
mean(temp)
|
|
temp <- numeric(nr)
|
|
set.seed(1234)
|
|
for (i_nr in seq_len(nr)) {
|
|
x <- rexp(n, rate = lambda)
|
|
temp_int <- conf_int(x)
|
|
temp[i_nr] <- ((temp_int[1] < 1 / lambda) & (1 / lambda < temp_int[2]))
|
|
}
|
|
mean(temp)
|
|
n <- 50
|
|
cat(paste("n = ", n, sep = ""), "\n")
|
|
temp <- numeric(nr)
|
|
set.seed(1234)
|
|
for (i_nr in seq_len(nr)) {
|
|
x <- rnorm(n, mu, sigma)
|
|
temp_int <- conf_int(x)
|
|
temp[i_nr] <- ((temp_int[1] < mu) & (mu < temp_int[2]))
|
|
}
|
|
mean(temp)
|
|
temp <- numeric(nr)
|
|
set.seed(1234)
|
|
for (i_nr in seq_len(nr)) {
|
|
x <- rchisq(n, df)
|
|
temp_int <- conf_int(x)
|
|
temp[i_nr] <- ((temp_int[1] < df) & (df < temp_int[2]))
|
|
}
|
|
mean(temp)
|
|
temp <- numeric(nr)
|
|
set.seed(1234)
|
|
for (i_nr in seq_len(nr)) {
|
|
x <- rexp(n, rate = lambda)
|
|
temp_int <- conf_int(x)
|
|
temp[i_nr] <- ((temp_int[1] < 1 / lambda) & (1 / lambda < temp_int[2]))
|
|
}
|
|
mean(temp)
|
|
n <- 100
|
|
cat(paste("n = ", n, sep = ""), "\n")
|
|
temp <- numeric(nr)
|
|
set.seed(1234)
|
|
for (i_nr in seq_len(nr)) {
|
|
x <- rnorm(n, mu, sigma)
|
|
temp_int <- conf_int(x)
|
|
temp[i_nr] <- ((temp_int[1] < mu) & (mu < temp_int[2]))
|
|
}
|
|
mean(temp)
|
|
temp <- numeric(nr)
|
|
set.seed(1234)
|
|
for (i_nr in seq_len(nr)) {
|
|
x <- rchisq(n, df)
|
|
temp_int <- conf_int(x)
|
|
temp[i_nr] <- ((temp_int[1] < df) & (df < temp_int[2]))
|
|
}
|
|
mean(temp)
|
|
temp <- numeric(nr)
|
|
set.seed(1234)
|
|
for (i_nr in seq_len(nr)) {
|
|
x <- rexp(n, rate = lambda)
|
|
temp_int <- conf_int(x)
|
|
temp[i_nr] <- ((temp_int[1] < 1 / lambda) & (1 / lambda < temp_int[2]))
|
|
}
|
|
mean(temp)
|