Data 604 Week 5 Discussion
Judd Anderman
September 27, 2017
Random number generation and testing
library(ggplot2)## Warning: package 'ggplot2' was built under R version 3.3.2library(gridExtra)lcg <- function(X_0, a, c, m, length) {
stream <- vector("numeric", length)
stream[1] <- X_0
for (i in 1:(length)) {
stream[i + 1] <- (a * stream[i] + c) %% m
}
stream <- stream[-1]
stream <- stream / m
return(stream)
}x1 <- lcg(27, 17, 43, 100, 1e4)
x2 <- lcg(123457, 7^5, 0, 2^31 - 1, 1e4)
pdf1 <- ggplot(as.data.frame(x1), aes(x1)) +
geom_line(stat = "density") +
coord_cartesian(xlim = c(0, 1)) +
labs(x = "x", y = "density")
ecdf1 <- ggplot(as.data.frame(x1), aes(x1)) +
stat_ecdf(geom = "step") +
coord_cartesian(xlim = c(0, 1)) +
geom_abline(slope = 1, intercept = 0,
color = "red", linetype = "dotted") +
labs(x = "x", y = "cumulative probability")
grid.arrange(pdf1, ecdf1, nrow = 1)
pdf2 <- ggplot(as.data.frame(x2), aes(x2)) +
geom_line(stat = "density") +
coord_cartesian(xlim = c(0, 1)) +
labs(x = "x", y = "density")
ecdf2 <- ggplot(as.data.frame(x2), aes(x2)) +
stat_ecdf(geom = "step") +
coord_cartesian(xlim = c(0, 1)) +
geom_abline(slope = 1, intercept = 0,
color = "red", linetype = "dotted") +
labs(x = "x", y = "cumulative probability")
grid.arrange(pdf2, ecdf2, nrow = 1)
df1 <- data.frame(x = x1)
df1$S <- sapply(df1$x, function(r) length(df1$x[df1$x <= r]) / length(df1$x))
head(df1)## x S
## 1 0.02 0.25
## 2 0.77 1.00
## 3 0.52 0.75
## 4 0.27 0.50
## 5 0.02 0.25
## 6 0.77 1.00df1 <- df1[order(df1$x), ]
rownames(df1) <- NULL
D_pos1 <- max(as.integer(rownames(df1)) / nrow(df1) - df1$x)
D_neg1 <- max(df1$x - ((as.integer(rownames(df1)) - 1) / nrow(df1)))
D1 <- max(D_pos1, D_neg1)
df2 <- data.frame(x = x2)
df2$S <- sapply(df2$x, function(r) length(df2$x[df2$x <= r]) / length(df2$x))
head(df2)## x S
## 1 0.96622007 0.9644
## 2 0.26071079 0.2596
## 3 0.76626223 0.7658
## 4 0.56933687 0.5662
## 5 0.84482919 0.8431
## 6 0.04426651 0.0425df2 <- df2[order(df2$x), ]
rownames(df2) <- NULL
D_pos2 <- max(as.integer(rownames(df2)) / nrow(df2) - df2$x)
D_neg2 <- max(df2$x - ((as.integer(rownames(df2)) - 1) / nrow(df2)))
D2 <- max(D_pos2, D_neg2)
ks.test(x1, "punif")## Warning in ks.test(x1, "punif"): ties should not be present for the
## Kolmogorov-Smirnov test##
## One-sample Kolmogorov-Smirnov test
##
## data: x1
## D = 0.23, p-value < 2.2e-16
## alternative hypothesis: two-sided# empirical cumulative distribution function is significantly different from continuous uniform cdf
ks.test(x2, "punif")##
## One-sample Kolmogorov-Smirnov test
##
## data: x2
## D = 0.0056387, p-value = 0.9082
## alternative hypothesis: two-sided# empirical cumulative distribution function is NOT significantly different from continuous uniform cdf
M1 <- floor(((length(x1) - 4) / 4) - 1)
corsum1 <- 0
for (k in 0:M1) {
corsum1 <- corsum1 + (x1[1 + k * 4] * x1[1 + (k + 1) * 4])
}
rho1 <- (1 / (M1 + 1)) * corsum1 - 0.25
sigma1 <- sqrt(13 * M1 + 7) / (12 * (M1 + 1))
(z1 <- rho1 / sigma1)## [1] -41.53148(half_alpha <- qnorm(0.975))## [1] 1.959964-1 * half_alpha <= z1 & z1 <= half_alpha## [1] FALSE# reject H_0 (i.e., autocorrelation with lag of 4 beginning at first random number = 0)
M2 <- floor(((length(x2) - 4) / 4) - 1)
corsum2 <- 0
for (k in 0:M2) {
corsum2 <- corsum2 + (x2[1 + k * 4] * x2[1 + (k + 1) * 4])
}
rho2 <- (1 / (M2 + 1)) * corsum2 - 0.25
sigma2 <- sqrt(13 * M2 + 7) / (12 * (M2 + 1))
(z2 <- rho2 / sigma2)## [1] 0.4691198(half_alpha <- qnorm(0.975))## [1] 1.959964-1 * half_alpha <= z2 & z2 <= half_alpha## [1] TRUE# fail to reject H_0
(ac1 <- acf(x1, type = "correlation", plot = FALSE))##
## Autocorrelations of series 'x1', by lag
##
## 0 1 2 3 4 5 6 7 8 9
## 1.000 -0.200 -0.600 -0.200 1.000 -0.200 -0.600 -0.200 0.999 -0.200
## 10 11 12 13 14 15 16 17 18 19
## -0.599 -0.200 0.999 -0.200 -0.599 -0.200 0.998 -0.200 -0.599 -0.200
## 20 21 22 23 24 25 26 27 28 29
## 0.998 -0.200 -0.599 -0.199 0.998 -0.200 -0.598 -0.199 0.997 -0.199
## 30 31 32 33 34 35 36 37 38 39
## -0.598 -0.199 0.997 -0.199 -0.598 -0.199 0.996 -0.199 -0.598 -0.199
## 40
## 0.996(ac2 <- acf(x2, type = "correlation", plot = FALSE))##
## Autocorrelations of series 'x2', by lag
##
## 0 1 2 3 4 5 6 7 8 9
## 1.000 -0.006 -0.022 -0.008 0.010 -0.016 -0.013 0.008 -0.003 -0.011
## 10 11 12 13 14 15 16 17 18 19
## -0.003 0.004 -0.001 -0.011 0.012 0.012 0.002 -0.001 -0.008 -0.020
## 20 21 22 23 24 25 26 27 28 29
## 0.025 -0.012 0.003 -0.005 0.015 0.023 0.018 0.000 0.007 -0.011
## 30 31 32 33 34 35 36 37 38 39
## -0.015 0.001 0.004 -0.018 -0.001 0.012 -0.010 0.002 0.011 0.016
## 40
## -0.003ggplot(data.frame(lag = ac1$lag, acf = ac1$acf), aes(lag, acf)) +
geom_point() +
geom_segment(aes(x = lag, xend = lag, y = 0, yend = acf)) +
geom_hline(yintercept = 0, linetype = "dashed") +
coord_cartesian(ylim = c(-1, 1))
ggplot(data.frame(lag = ac2$lag, acf = ac2$acf), aes(lag, acf)) +
geom_point() +
geom_segment(aes(x = lag, xend = lag, y = 0, yend = acf)) +
geom_hline(yintercept = 0, linetype = "dashed") +
coord_cartesian(ylim = c(-1, 1))
Kelton et al. (2014): Problem 4.13
lambda <- 100 / 60
mu <- 50 / 60
c <- 3
mmc <- function(lambda, mu, c) {
p <- lambda / (c * mu)
f <- function (c, p) {
x <- 0
for (n in 0:(c - 1)) {
x <- x + ((c * p)^n / factorial(n))
}
return(x)
}
p_0 <- 1 / ((c * p)^c / (factorial(c) * (1 - p)) + f(c, p))
L_q <- (p * (c * p)^c * p_0) / (factorial(c) * (1 - p)^2)
W_q <- L_q / lambda
L <- L_q + (lambda / mu)
W <- L / lambda
return(data.frame(W_q, W, L_q, L, p))
}
(MM3 <- mmc(lambda, mu, c))## W_q W L_q L p
## 1 0.5333333 1.733333 0.8888889 2.888889 0.6666667
Simio model
300 hour single run simulation results