## ----setup, include = FALSE--------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ## ----fig.width = 6, fig.height = 4, fig.show = 'hold', echo = FALSE, warning = FALSE, results = "hide"---- library(conf) nukedata <- c(1728, 1986, 10746) x1 <- crplot(nukedata, 0.1, "gamma", sf = c(2, 2), info = TRUE, repair = FALSE, showplot = FALSE) x2 <- crplot(nukedata, 0.1, "gamma", sf = c(2, 2), jumpinfo = TRUE, repair = TRUE, showplot = FALSE) jumpxy <- x2$repair$q4jumpxy gapindex <- which(x1$theta == max(x1$theta)) gapl <- c(x1$theta[gapindex + 1], x1$kappa[gapindex + 1]) gapr <- c(x1$theta[gapindex ], x1$kappa[gapindex ]) par(mar = c(4, 4, 1, 2), xpd = TRUE) plot(c(x2$theta, x2$theta[1]), c(x2$kappa, x2$kappa[1]), xlab = expression(theta), ylab = expression(kappa), type = "l", axes = FALSE, col = "white") label1 <- format(c(min(x2$theta), x2$thetahat, max(x1$theta), max(x2$theta)), digits = 3) label2 <- format(c(min(x2$kappa), x2$kappahat, max(x2$kappa)), digits = 3) axis(side = 1, at = label1, labels = TRUE) axis(side = 2, at = label2, labels = TRUE, las = 1) # identify points created from jumpxy jump-center newpoints <- intersect(which(x2$theta > gapl[1]), which (x2$kappa > gapr[2])) # draw angles taken from jump center to ID new points m2 <- (x2$kappa[newpoints] - jumpxy[2]) / (x2$theta[newpoints] - jumpxy[1]) par(xpd = FALSE) segments(jumpxy[1], jumpxy[2], jumpxy[1] + (max(x2$kappa) - jumpxy[2]) / m2, max(x2$kappa), col = "gray75", lwd = 0.5) par(xpd = TRUE) points(x2$theta[newpoints], x2$kappa[newpoints], pch = 1) points(x2$theta[newpoints], x2$kappa[newpoints], pch = 16, cex = 0.5, col = "blue") # draw CR boundary and its "new" points lines(c(x2$theta, x2$theta[1]), c(x2$kappa, x2$kappa[1])) points(x2$thetahat, x2$kappahat, pch = 3) # identify the last CR point preceeding "gap" in its perimeter # and draw a line from the MLE to this point gapindex <- which(x1$theta == max(x1$theta)) gapl <- c(x1$theta[gapindex + 1], x1$kappa[gapindex + 1]) gapr <- c(x1$theta[gapindex ], x1$kappa[gapindex ]) # show dotted line from MLE along edge of CR inaccessible area m <- (x1$kappahat - gapl[2]) / (x1$thetahat - gapl[1]) # slope of line marking inaccessible region lines(c(x1$thetahat + 1600, gapr[1] + 8000), c(x1$kappahat + 1600 * m, gapr[2] + 8000 * m), lty = 3, col = "red", lwd = 1.5) # make additional markings halving distance below the left-edge of the gap half <- (x1$kappa[gapindex + 1] - x1$kappa[gapindex]) / 2 xcol <- "gray30" #lines(c(x1$theta[gapindex + 1] - 700, x1$theta[gapindex + 1] + 700), min(x1$kappa) + c(half, half), col = xcol) arrows(gapl[1], gapl[2] - 0.1, gapl[1], gapl[2] - half + 0.005, code = 3, length = 0.04, col = "blue", lwd = 1) text(gapl[1] + 5200, gapl[2] - 0.32, "0.5(gap)", col = "blue", lwd = 1, cex = 0.8) # show arrow to new jump center, and that center point # obtained jump center location (jumpx, jumpy) by unique, targeted, conf modifications m2 <- (x1$kappahat - jumpxy[2]) / (x1$thetahat - jumpxy[1]) arrows(x1$thetahat + 1600, x1$kappahat + 1600 * m2, jumpxy[1] - 2000, jumpxy[2] - 2000 * m2, length = 0.1, lty = 1, angle = 15, lwd = 0.9, col = "blue") points(jumpxy[1], jumpxy[2], pch = 24, bg = "red") text(x = x1$theta[gapindex + 1] - 1000, y = 0.6, '{', srt = 0, cex = 2.4, family = 'Helvetica Neue UltraLight', col = "purple", lwd = 1.5) text(x = x1$theta[gapindex + 1] - 4500, y = 0.55, 'gap', col = "purple", cex = 0.9) points(c(gapl[1], gapr[1]), c(gapl[2], gapr[2]), pch = 24, bg = c("green", "yellow")) legend("topright", legend = c("jump-center (JC)", "JC left boundary", "JC right boundary", "JC repair points"), pch = c(rep(24, 3), 21), pt.bg = c("red", "green", "yellow", "blue"), pt.cex = c(1, 1, 1, 0.7), bty = "n", cex = 0.8) ## ----------------------------------------------------------------------------- bc_obs <- c(230, 334, 423, 990, 1009, 1510) bc_cen <- c(rep(50, 288), rep(150, 148), rep(250, 124), rep(350, 111), rep(450, 106), rep(550, 99), rep(650, 110), rep(750, 114), rep(850, 119), rep(950, 127), rep(1050, 123), rep(1150, 93), rep(1250, 47), rep(1350, 41), rep(1450, 27), rep(1550, 11), rep(1650, 6), rep(1850, 1), rep(2050, 2)) bc <- c(bc_obs, bc_cen) cen <- c(rep(1, length(bc_obs)), rep(0, length(bc_cen))) print(length(bc)) ## ----fig.width = 7, fig.height = 5, fig.show = 'hold', warning = FALSE, eval = FALSE---- # distns <- c("cauchy", "gamma", "llogis", "logis", "norm", "weibull") # par(mfrow = c(2, 3)) # plot in a 2-row, 3-column grid # for (i in 1:length(distns)) { # try(crplot(dataset = bc, alpha = 0.1, distn = distns[i], cen = cen, # sf = c(0, 2), ylas = 1, main = distns[i])) # } ## ----fig.width = 7, fig.height = 5, fig.show = 'hold', warning = FALSE, echo = FALSE---- # results of above chunk are input directly here to save processing time on vignette (otherwise prohibitively long for CRAN) distns <- c("cauchy", "gamma", "llogis", "logis", "norm", "weibull") par(mfrow = c(2, 3)) # plot in a 2-row, 3-column grid for (i in 1:length(distns)) { if (i != 2) { try(crplot(dataset = bc, alpha = 0.1, distn = distns[i], cen = cen, sf = c(0, 2), ylas = 1, main = distns[i])) } else if (i == 2) { a <- NULL a$theta1 <- c(11092.63, 10433.13, 9807.861, 8634.626, 8073.472, 7796.363, 7520.319, 6952.024, 6371.019, 5763.591, 5061.823, 4510.083, 4126.239, 4077.722, 4069.232, 4065.99, 4065.685, 4065.64, 4065.629, 4065.625, 3988.346, 3894.447, 3810.225, 3733.788, 3663.776, 3599.174, 3539.203, 3483.248, 3430.814, 3381.497, 3334.963, 3290.929, 3249.157, 3209.441, 3171.605, 3135.494, 3100.973, 3067.922, 2422.232, 2142.176, 1874.156, 1867.841, 1861.526, 1855.118, 1848.709, 1842.203, 1835.697, 1829.092, 1822.485, 1815.777, 1809.066, 1802.251, 1795.434, 1788.509, 1781.58, 1774.541, 1767.499, 1760.342, 1753.18, 1745.901, 1738.617, 1731.212, 1723.8, 1716.264, 1708.721, 1701.048, 1693.369, 1685.556, 1677.735, 1669.776, 1661.809, 1653.699, 1645.58, 1637.314, 1629.039, 1620.611, 1612.174, 1603.581, 1594.978, 1586.215, 1577.441, 1568.504, 1559.558, 1550.446, 1541.325, 1532.037, 1522.744, 1513.282, 1503.82, 1494.193, 1484.572, 1474.792, 1465.028, 1455.118, 1445.239, 1425.278, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.423, 1557.423, 1557.424, 1557.425, 1557.426, 1557.429, 1557.433, 1557.437, 1557.444, 1557.45, 1557.469, 1557.489, 1557.527, 1557.603, 1557.723, 1557.844, 1558.102, 1558.36, 1558.877, 1559.69, 1560.503, 1561.782, 1563.064, 1565.083, 1567.108, 1578.992, 1591.141, 1616.27, 1653.549, 1693.321, 1753.816, 1820.462, 2097.608, 2518.007, 2826.445, 3253.11, 3937.509, 4439.756, 5098.404, 5620.529, 6212.872, 6854.903, 7520.319, 8305.757, 9067.957, 9792.367, 10483.48, 11730.83, 12880.77, 14784.04, 16510.45, 19141.57, 21541.62, 25392.41, 27192.46, 28930.35, 35822.83, 42440.97, 48878.38, 49595.8, 50309.54, 51411.59, 51837.43, 52262.06, 52591.09, 52919.42, 53173.97, 53428.11, 53625.22, 53822.09, 53974.84, 54127.44, 54245.87, 54364.21, 54456.07, 54547.88, 54619.15, 54690.39, 54745.7, 54801, 54843.93, 54865.39, 54886.85, 54894.62, 54902.38, 54908.34, 54914.29, 54918.86, 54923.43, 54926.93, 54930.43, 54933.12, 54935.81, 54937.87, 54939.93, 54941.51, 54943.09, 54944.31, 54945.52, 339933.9, 343422.4, 346896.7, 350350.8, 353777.2, 355481.7, 356326.7, 356747.3, 357166.7, 357542.7, 357917.7, 358430, 181819.7, 169469.9, 158491.6, 150581.6, 144265.7, 138960.9, 134366.9, 130305.5, 126660.7, 123352.2, 120322.1, 117526.7, 114932.2, 112512, 110244.6, 108112.4, 106100.7, 104197.1, 102391.4, 100674.4, 99038.36, 97476.61, 95983.17, 94552.8, 93180.87, 82254.1, 72703.74, 68327.03, 64176.95, 56455.64, 52887.15, 49455.84, 46147.03, 42944.25, 39851.51, 36832.26, 33872.02, 30943.79, 30708.99, 30474.16, 29996.34, 29509.8, 29013.94, 28508.08, 27991.44, 27463.14, 26922.14, 26367.22, 25796.96, 25209.63, 24603.14, 23974.9, 23321.66, 22639.19, 21951.15, 21922.4, 21223.25, 18865.36, 17491.4, 16398.69, 15484.77, 13878.77, 12551.42) a$theta2 <- c(2.142689, 2.180192, 2.218991, 2.302236, 2.347989, 2.372294, 2.397765, 2.454672, 2.520183, 2.598519, 2.705241, 2.805313, 2.885938, 2.896894, 2.898831, 2.899572, 2.899642, 2.899652, 2.899655, 2.899655, 2.917575, 2.940029, 2.960838, 2.9803, 2.998632, 3.015996, 3.032517, 3.048296, 3.063412, 3.077932, 3.091912, 3.105397, 3.118428, 3.131039, 3.14326, 3.155117, 3.166634, 3.177831, 3.437111, 3.580792, 3.742568, 3.746705, 3.750857, 3.755087, 3.759334, 3.763661, 3.768005, 3.772433, 3.776879, 3.781411, 3.785962, 3.790602, 3.795263, 3.800016, 3.80479, 3.809659, 3.814551, 3.819543, 3.824557, 3.829675, 3.834816, 3.840065, 3.845339, 3.850723, 3.856134, 3.861659, 3.867212, 3.872884, 3.878583, 3.884406, 3.890258, 3.896237, 3.902244, 3.908384, 3.914552, 3.920856, 3.927188, 3.933659, 3.940158, 3.946797, 3.953463, 3.960271, 3.967102, 3.974075, 3.981066, 3.988197, 3.995339, 4.002615, 4.009891, 4.017291, 4.024677, 4.032168, 4.039625, 4.047163, 4.054636, 4.069564, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502057, 3.502056, 3.502055, 3.502054, 3.502052, 3.502049, 3.502043, 3.502038, 3.502029, 3.502021, 3.501996, 3.501971, 3.501922, 3.501822, 3.501666, 3.501509, 3.501173, 3.500838, 3.500166, 3.499112, 3.498058, 3.496403, 3.494747, 3.492148, 3.489547, 3.474448, 3.459287, 3.428771, 3.385445, 3.341559, 3.278935, 3.215022, 2.993232, 2.747052, 2.609885, 2.459118, 2.278387, 2.176874, 2.069931, 2.000331, 1.933303, 1.871587, 1.816841, 1.761516, 1.715315, 1.676807, 1.644061, 1.592815, 1.552607, 1.497032, 1.455406, 1.40337, 1.364588, 1.314333, 1.294595, 1.277305, 1.221507, 1.181113, 1.149826, 1.146713, 1.14368, 1.139119, 1.137395, 1.135696, 1.134393, 1.133105, 1.132114, 1.131132, 1.130375, 1.129623, 1.129042, 1.128464, 1.128018, 1.127573, 1.127228, 1.126885, 1.126619, 1.126353, 1.126147, 1.125942, 1.125783, 1.125703, 1.125624, 1.125595, 1.125566, 1.125544, 1.125522, 1.125505, 1.125488, 1.125475, 1.125462, 1.125453, 1.125443, 1.125435, 1.125427, 1.125421, 1.125416, 1.125411, 1.125407, 0.8874298, 0.8872036, 0.8870392, 0.8869431, 0.886923, 0.8869445, 0.8869636, 0.8869753, 0.8869885, 0.8870016, 0.887016, 0.8870376, 1.090469, 1.107421, 1.123776, 1.136433, 1.147137, 1.156582, 1.165129, 1.172993, 1.180315, 1.187192, 1.193698, 1.199885, 1.205795, 1.211463, 1.216914, 1.222172, 1.227255, 1.232179, 1.236957, 1.241602, 1.246122, 1.250528, 1.254826, 1.259024, 1.263128, 1.298867, 1.335565, 1.354559, 1.374107, 1.415339, 1.437012, 1.459772, 1.483813, 1.509392, 1.536653, 1.56617, 1.598478, 1.634452, 1.637537, 1.640654, 1.6471, 1.653809, 1.660804, 1.66811, 1.675758, 1.683781, 1.69222, 1.70112, 1.710539, 1.720543, 1.731216, 1.742661, 1.755008, 1.768427, 1.782534, 1.783137, 1.79814, 1.854264, 1.891683, 1.924498, 1.954384, 2.013396, 2.069931) a$theta1hat <- 7520.318946 a$theta2hat <- 2.069931 plot(c(a$theta1, a$theta1[1]), c(a$theta2, a$theta2[1]), type = "l", axes = FALSE, xlab = expression(theta), ylab = expression(kappa), main = "gamma") points(a$theta1, a$theta2, lwd = 0.65) points(a$theta1hat, a$theta2hat, pch = 3) axis(side = 1, at = format(c(range(a$theta1), a$theta1hat), digits = 2)) axis(side = 2, at = format(c(range(a$theta2), a$theta2hat), digits = 2), las = 1) } } ## ----fig.show = 'hold', warning = FALSE, eval = FALSE------------------------- # x <- crplot(bc, alpha = 0.1, distn = "gamma", cen = cen, jumpinfo = TRUE, showjump = TRUE, # main = "with repairs") # #> [1] "Confidence region plot complete; made using 276 boundary points." # crplot(bc, alpha = 0.1, distn = "gamma", cen = cen, repair = FALSE, # main = "without repairs") # #> [1] "Confidence region plot complete; made using 134 boundary points." # jumppoints <- t(matrix(c(x$repair$q2jumpxy, x$repair$q2jumpL, x$repair$q2jumpR, # x$repair$q4jumpxy, x$repair$q4jumpL, x$repair$q4jumpR), nrow = 2)) # points(jumppoints, pch = 24, bg = rep(c("red", "green", "yellow"))) # print(labels(x$repair)) # #> [1] "q2jumpuphill" "q4jumpuphill" "q2jumpshift" "q4jumpshift" # #> [5] "q2jumpxy" "q4jumpxy" "q2jumpL" "q4jumpL" # #> [9] "q2jumpR" "q4jumpR" "q2gaptype" "q4gaptype" ## ----fig.show = 'hold', warning = FALSE, results='hide', echo = FALSE--------- # results of above chunk are input directly here to save processing time on vignette (otherwise prohibitively long for CRAN) par(mar = c(4, 4.5, 2, 1.5)) a$theta1 <- c(11092.63, 10433.13, 9807.861, 8634.626, 8073.472, 7796.363, 7520.319, 6952.024, 6371.019, 5763.591, 5061.823, 4510.083, 4126.239, 4077.722, 4069.232, 4065.99, 4065.685, 4065.64, 4065.629, 4065.625, 3988.346, 3894.447, 3810.225, 3733.788, 3663.776, 3599.174, 3539.203, 3483.248, 3430.814, 3381.497, 3334.963, 3290.929, 3249.157, 3209.441, 3171.605, 3135.494, 3100.973, 3067.922, 2422.232, 2142.176, 1874.156, 1867.841, 1861.526, 1855.118, 1848.709, 1842.203, 1835.697, 1829.092, 1822.485, 1815.777, 1809.066, 1802.251, 1795.434, 1788.509, 1781.58, 1774.541, 1767.499, 1760.342, 1753.18, 1745.901, 1738.617, 1731.212, 1723.8, 1716.264, 1708.721, 1701.048, 1693.369, 1685.556, 1677.735, 1669.776, 1661.809, 1653.699, 1645.58, 1637.314, 1629.039, 1620.611, 1612.174, 1603.581, 1594.978, 1586.215, 1577.441, 1568.504, 1559.558, 1550.446, 1541.325, 1532.037, 1522.744, 1513.282, 1503.82, 1494.193, 1484.572, 1474.792, 1465.028, 1455.118, 1445.239, 1425.278, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.423, 1557.423, 1557.424, 1557.425, 1557.426, 1557.429, 1557.433, 1557.437, 1557.444, 1557.45, 1557.469, 1557.489, 1557.527, 1557.603, 1557.723, 1557.844, 1558.102, 1558.36, 1558.877, 1559.69, 1560.503, 1561.782, 1563.064, 1565.083, 1567.108, 1578.992, 1591.141, 1616.27, 1653.549, 1693.321, 1753.816, 1820.462, 2097.608, 2518.007, 2826.445, 3253.11, 3937.509, 4439.756, 5098.404, 5620.529, 6212.872, 6854.903, 7520.319, 8305.757, 9067.957, 9792.367, 10483.48, 11730.83, 12880.77, 14784.04, 16510.45, 19141.57, 21541.62, 25392.41, 27192.46, 28930.35, 35822.83, 42440.97, 48878.38, 49595.8, 50309.54, 51411.59, 51837.43, 52262.06, 52591.09, 52919.42, 53173.97, 53428.11, 53625.22, 53822.09, 53974.84, 54127.44, 54245.87, 54364.21, 54456.07, 54547.88, 54619.15, 54690.39, 54745.7, 54801, 54843.93, 54865.39, 54886.85, 54894.62, 54902.38, 54908.34, 54914.29, 54918.86, 54923.43, 54926.93, 54930.43, 54933.12, 54935.81, 54937.87, 54939.93, 54941.51, 54943.09, 54944.31, 54945.52, 339933.9, 343422.4, 346896.7, 350350.8, 353777.2, 355481.7, 356326.7, 356747.3, 357166.7, 357542.7, 357917.7, 358430, 181819.7, 169469.9, 158491.6, 150581.6, 144265.7, 138960.9, 134366.9, 130305.5, 126660.7, 123352.2, 120322.1, 117526.7, 114932.2, 112512, 110244.6, 108112.4, 106100.7, 104197.1, 102391.4, 100674.4, 99038.36, 97476.61, 95983.17, 94552.8, 93180.87, 82254.1, 72703.74, 68327.03, 64176.95, 56455.64, 52887.15, 49455.84, 46147.03, 42944.25, 39851.51, 36832.26, 33872.02, 30943.79, 30708.99, 30474.16, 29996.34, 29509.8, 29013.94, 28508.08, 27991.44, 27463.14, 26922.14, 26367.22, 25796.96, 25209.63, 24603.14, 23974.9, 23321.66, 22639.19, 21951.15, 21922.4, 21223.25, 18865.36, 17491.4, 16398.69, 15484.77, 13878.77, 12551.42) a$theta2 <- c(2.142689, 2.180192, 2.218991, 2.302236, 2.347989, 2.372294, 2.397765, 2.454672, 2.520183, 2.598519, 2.705241, 2.805313, 2.885938, 2.896894, 2.898831, 2.899572, 2.899642, 2.899652, 2.899655, 2.899655, 2.917575, 2.940029, 2.960838, 2.9803, 2.998632, 3.015996, 3.032517, 3.048296, 3.063412, 3.077932, 3.091912, 3.105397, 3.118428, 3.131039, 3.14326, 3.155117, 3.166634, 3.177831, 3.437111, 3.580792, 3.742568, 3.746705, 3.750857, 3.755087, 3.759334, 3.763661, 3.768005, 3.772433, 3.776879, 3.781411, 3.785962, 3.790602, 3.795263, 3.800016, 3.80479, 3.809659, 3.814551, 3.819543, 3.824557, 3.829675, 3.834816, 3.840065, 3.845339, 3.850723, 3.856134, 3.861659, 3.867212, 3.872884, 3.878583, 3.884406, 3.890258, 3.896237, 3.902244, 3.908384, 3.914552, 3.920856, 3.927188, 3.933659, 3.940158, 3.946797, 3.953463, 3.960271, 3.967102, 3.974075, 3.981066, 3.988197, 3.995339, 4.002615, 4.009891, 4.017291, 4.024677, 4.032168, 4.039625, 4.047163, 4.054636, 4.069564, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502057, 3.502056, 3.502055, 3.502054, 3.502052, 3.502049, 3.502043, 3.502038, 3.502029, 3.502021, 3.501996, 3.501971, 3.501922, 3.501822, 3.501666, 3.501509, 3.501173, 3.500838, 3.500166, 3.499112, 3.498058, 3.496403, 3.494747, 3.492148, 3.489547, 3.474448, 3.459287, 3.428771, 3.385445, 3.341559, 3.278935, 3.215022, 2.993232, 2.747052, 2.609885, 2.459118, 2.278387, 2.176874, 2.069931, 2.000331, 1.933303, 1.871587, 1.816841, 1.761516, 1.715315, 1.676807, 1.644061, 1.592815, 1.552607, 1.497032, 1.455406, 1.40337, 1.364588, 1.314333, 1.294595, 1.277305, 1.221507, 1.181113, 1.149826, 1.146713, 1.14368, 1.139119, 1.137395, 1.135696, 1.134393, 1.133105, 1.132114, 1.131132, 1.130375, 1.129623, 1.129042, 1.128464, 1.128018, 1.127573, 1.127228, 1.126885, 1.126619, 1.126353, 1.126147, 1.125942, 1.125783, 1.125703, 1.125624, 1.125595, 1.125566, 1.125544, 1.125522, 1.125505, 1.125488, 1.125475, 1.125462, 1.125453, 1.125443, 1.125435, 1.125427, 1.125421, 1.125416, 1.125411, 1.125407, 0.8874298, 0.8872036, 0.8870392, 0.8869431, 0.886923, 0.8869445, 0.8869636, 0.8869753, 0.8869885, 0.8870016, 0.887016, 0.8870376, 1.090469, 1.107421, 1.123776, 1.136433, 1.147137, 1.156582, 1.165129, 1.172993, 1.180315, 1.187192, 1.193698, 1.199885, 1.205795, 1.211463, 1.216914, 1.222172, 1.227255, 1.232179, 1.236957, 1.241602, 1.246122, 1.250528, 1.254826, 1.259024, 1.263128, 1.298867, 1.335565, 1.354559, 1.374107, 1.415339, 1.437012, 1.459772, 1.483813, 1.509392, 1.536653, 1.56617, 1.598478, 1.634452, 1.637537, 1.640654, 1.6471, 1.653809, 1.660804, 1.66811, 1.675758, 1.683781, 1.69222, 1.70112, 1.710539, 1.720543, 1.731216, 1.742661, 1.755008, 1.768427, 1.782534, 1.783137, 1.79814, 1.854264, 1.891683, 1.924498, 1.954384, 2.013396, 2.069931) a$theta1hat <- 7520.319 # MLE a$theta2hat <- 2.069931 # MLE pJ <- matrix(c(1874.156, 3.177831, 30943.79, 1.263128), nrow = 2, byrow = TRUE) # jump-center location pL <- matrix(c(1557.422, 3.502058, 21922.4, 1.783137), nrow = 2, byrow = TRUE) # jump-center left inaccessible region boundary pR <- matrix(c(4065.629, 2.899655, 54943.09, 1.125416), nrow = 2, byrow = TRUE) # jump-center right inaccessible region boundary #---> this point forward is universal for above inputs plot(c(a$theta1, a$theta1[1]), c(a$theta2, a$theta2[1]), type = "l", axes = FALSE, xlab = expression(theta), ylab = expression(kappa), main = "with repairs") # ensure axis labels are accurate points(a$theta1, a$theta2, lwd = 0.65) axis(side = 1, at = c(a$theta1hat, range(a$theta1))) axis(side = 2, at = c(a$theta2hat, range(a$theta2))) for (i in 1:nrow(pL)) { # ID & plot jump-center repair points in each quad m1R <- abs(a$theta1 - pR[i, 1]) / diff(range(a$theta1)) # match theta1 values m2R <- abs(a$theta2 - pR[i, 2]) / diff(range(a$theta2)) # match theta2 values m1L <- abs(a$theta1 - pL[i, 1]) / diff(range(a$theta1)) # match theta1 values m2L <- abs(a$theta2 - pL[i, 2]) / diff(range(a$theta2)) # match theta2 values fm <- which.min(m1R + m2R) + 1 # "right" points (nearest right boundary point) to <- which.min(m1L + m2L) - 1 # "to" points (nearest left boundary point) points(a$theta1[fm:to], a$theta2[fm:to], pch = 16, cex = 0.4, col = "blue") } points(pL, pch = 24, bg = "green") # points on left-side of inaccessible region points(pR, pch = 24, bg = "yellow") # points on right-side of inaccessible region points(pJ, pch = 24, bg = "red") # jump-center legend("topright", legend = c("CR boundary points", "jump-center (JC)", "JC left boundary", "JC right boundary", "JC repair points"), pch = c(1, rep(24, 3), 21), pt.bg = c("black", "red", "green", "yellow", "blue"), pt.cex = c(1, 1, 1, 1, 0.7), bty = "n", cex = 0.8) crplot(bc, alpha = 0.1, distn = "gamma", cen = cen, repair = FALSE, main = "without repairs") points(pJ, pch = 24, bg = "red") points(pL, pch = 24, bg = "green") points(pR, pch = 24, bg = "yellow") ## ----fig.width = 3.5, fig.height = 3.5, fig.show = 'hold', warning = FALSE, eval = FALSE---- # crplot(bc, alpha = 0.1, distn = "gamma", cen = cen, jumpuphill = 0.1, showjump = TRUE) # #> [1] "Confidence region plot complete; made using 277 boundary points." ## ----fig.width = 3.5, fig.height = 3.5, fig.show = 'hold', echo = FALSE------- par(mar = c(4, 4.5, 2, 1.5)) # results of above chunk are input directly here to save processing time on vignette (otherwise prohibitively long for CRAN) a$theta1 <- c(11092.63, 10433.13, 9807.861, 8634.626, 8073.472, 7796.363, 7520.319, 6952.024, 6371.019, 5763.591, 5061.823, 4510.083, 4126.239, 4077.722, 4069.232, 4065.99, 4065.685, 4065.64, 4065.629, 4065.625, 3525.655, 3140.348, 2525.686, 2256.465, 1998.091, 1935.64, 1873.007, 1799.155, 1724.441, 1633.425, 1540.458, 1428.438, 1328.577, 1309.528, 1301.099, 1297.174, 1295.285, 1293.445, 1292.993, 1292.544, 1291.656, 1291.216, 1290.78, 1290.346, 1290.131, 1289.916, 1289.702, 1289.489, 1289.276, 1289.065, 1288.854, 1288.644, 1288.435, 1288.33, 1288.226, 1288.122, 1288.019, 1287.915, 1287.812, 1287.708, 1287.606, 1287.503, 1287.4, 1287.196, 1286.992, 1286.789, 1286.587, 1286.385, 1286.185, 1285.985, 1285.786, 1285.39, 1284.997, 1284.608, 1284.221, 1283.838, 1283.458, 1281.969, 1281.244, 1280.53, 1277.803, 1275.273, 1270.803, 1260.504, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.423, 1557.423, 1557.424, 1557.425, 1557.426, 1557.429, 1557.433, 1557.437, 1557.444, 1557.45, 1557.469, 1557.489, 1557.527, 1557.603, 1557.723, 1557.844, 1558.102, 1558.36, 1558.877, 1559.69, 1560.503, 1561.782, 1563.064, 1565.083, 1567.108, 1578.992, 1591.141, 1616.27, 1653.549, 1693.321, 1753.816, 1820.462, 2097.608, 2518.007, 2826.445, 3253.11, 3937.509, 4439.756, 5098.404, 5620.529, 6212.872, 6854.903, 7520.319, 8305.757, 9067.957, 9792.367, 10483.48, 11730.83, 12880.77, 14784.04, 16510.45, 19141.57, 21541.62, 25392.41, 27192.46, 28930.35, 35822.83, 42440.97, 48878.38, 49595.8, 50309.54, 51411.59, 51837.43, 52262.06, 52591.09, 52919.42, 53173.97, 53428.11, 53625.22, 53822.09, 53974.84, 54127.44, 54245.87, 54364.21, 54456.07, 54547.88, 54619.15, 54690.39, 54745.7, 54801, 54843.93, 54865.39, 54886.85, 54894.62, 54902.38, 54908.34, 54914.29, 54918.86, 54923.43, 54926.93, 54930.43, 54933.12, 54935.81, 54937.87, 54939.93, 54941.51, 54943.09, 54944.31, 54945.52, 58267.4, 67489.96, 77564.88, 88465.32, 100202.3, 112805.2, 126313, 140769.7, 156222.2, 172718.2, 190304.7, 209025.9, 228918.9, 250008.4, 272294.4, 284150.4, 289973.6, 295728.1, 315616.9, 323398.6, 331023.1, 332573.3, 334116.1, 335260.1, 335830.6, 336400, 336625.5, 336850.9, 337016.2, 337098.8, 337140.1, 337160.7, 337171, 337181.3, 337182, 337182.3, 337182.4, 337182.6, 337182.6, 337182.7, 337182.7, 337182.7, 337182.7, 337182.7, 337182.7, 337182.7, 337182.7, 337182.7, 337182.7, 337182.7, 337182.7, 142726.3, 142726.3, 142726.3, 142726.2, 142725.4, 142722.6, 142708.3, 142655.4, 142459, 141740.2, 137208.2, 116687.7, 93851.46, 81356.56, 72522.87, 64731.15, 61139.05, 57720.12, 51325.18, 45474.96, 39998.33, 37382.62, 34822.38, 29813.07, 27877.58, 25926.37, 23468.09, 22194.71, 21951.15, 21922.4, 21223.25, 18865.36, 17491.4, 16398.69, 15484.77, 13878.77, 12551.42) a$theta2 <- c(2.142689, 2.180192, 2.218991, 2.302236, 2.347989, 2.372294, 2.397765, 2.454672, 2.520183, 2.598519, 2.705241, 2.805313, 2.885938, 2.896894, 2.898831, 2.899572, 2.899642, 2.899652, 2.899655, 2.899655, 3.036305, 3.153512, 3.389528, 3.519332, 3.664494, 3.703103, 3.74332, 3.792716, 3.844882, 3.91128, 3.981732, 4.067219, 4.133202, 4.141432, 4.14397, 4.144827, 4.145151, 4.145405, 4.145458, 4.145506, 4.145589, 4.145624, 4.145655, 4.145681, 4.145693, 4.145703, 4.145713, 4.145721, 4.145728, 4.145734, 4.145739, 4.145743, 4.145746, 4.145747, 4.145747, 4.145748, 4.145748, 4.145748, 4.145747, 4.145747, 4.145746, 4.145745, 4.145743, 4.145739, 4.145734, 4.145728, 4.145721, 4.145713, 4.145704, 4.145693, 4.145682, 4.145656, 4.145625, 4.14559, 4.145551, 4.145507, 4.145459, 4.145223, 4.14508, 4.144918, 4.144101, 4.143011, 4.140027, 4.082455, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502057, 3.502056, 3.502055, 3.502054, 3.502052, 3.502049, 3.502043, 3.502038, 3.502029, 3.502021, 3.501996, 3.501971, 3.501922, 3.501822, 3.501666, 3.501509, 3.501173, 3.500838, 3.500166, 3.499112, 3.498058, 3.496403, 3.494747, 3.492148, 3.489547, 3.474448, 3.459287, 3.428771, 3.385445, 3.341559, 3.278935, 3.215022, 2.993232, 2.747052, 2.609885, 2.459118, 2.278387, 2.176874, 2.069931, 2.000331, 1.933303, 1.871587, 1.816841, 1.761516, 1.715315, 1.676807, 1.644061, 1.592815, 1.552607, 1.497032, 1.455406, 1.40337, 1.364588, 1.314333, 1.294595, 1.277305, 1.221507, 1.181113, 1.149826, 1.146713, 1.14368, 1.139119, 1.137395, 1.135696, 1.134393, 1.133105, 1.132114, 1.131132, 1.130375, 1.129623, 1.129042, 1.128464, 1.128018, 1.127573, 1.127228, 1.126885, 1.126619, 1.126353, 1.126147, 1.125942, 1.125783, 1.125703, 1.125624, 1.125595, 1.125566, 1.125544, 1.125522, 1.125505, 1.125488, 1.125475, 1.125462, 1.125453, 1.125443, 1.125435, 1.125427, 1.125421, 1.125416, 1.125411, 1.125407, 1.113641, 1.085536, 1.060603, 1.038452, 1.018668, 1.000895, 0.9848447, 0.9702901, 0.9570539, 0.9449997, 0.9340263, 0.9240639, 0.9150736, 0.9070519, 0.9000417, 0.8968784, 0.8954648, 0.8941586, 0.8903587, 0.8891941, 0.8882506, 0.8880848, 0.8879289, 0.8878195, 0.887767, 0.8877159, 0.887696, 0.8876764, 0.8876621, 0.887655, 0.8876515, 0.8876497, 0.8876488, 0.887648, 0.8876479, 0.8876479, 0.8876479, 0.8876478, 0.8876478, 0.8876478, 0.8876478, 0.8876478, 0.8876478, 0.8876478, 0.8876478, 0.8876478, 0.8876478, 0.8876478, 0.8876478, 0.8876478, 0.8876478, 1.149833, 1.149833, 1.149833, 1.149834, 1.149835, 1.14984, 1.149865, 1.149959, 1.150305, 1.151579, 1.159801, 1.201778, 1.261112, 1.302074, 1.33632, 1.371402, 1.389505, 1.408091, 1.447121, 1.488979, 1.535295, 1.560551, 1.587701, 1.649609, 1.67747, 1.708376, 1.752199, 1.777472, 1.782534, 1.783137, 1.79814, 1.854264, 1.891683, 1.924498, 1.954384, 2.013396, 2.069931) a$theta1hat <- 7520.319 # MLE a$theta2hat <- 2.069931 # MLE pJ <- matrix(c(1998.091, 3.153512, 29813.07, 1.302074), nrow = 2, byrow = TRUE) # jump-center location pL <- matrix(c(1557.422, 3.502058, 21922.4, 1.783137), nrow = 2, byrow = TRUE) # jump-center left inaccessible region boundary pR <- matrix(c(4065.629, 2.899655, 54943.09, 1.125416), nrow = 2, byrow = TRUE) # jump-center right inaccessible region boundary #---> this point forward is universal for above inputs plot(c(a$theta1, a$theta1[1]), c(a$theta2, a$theta2[1]), type = "l", axes = FALSE, xlab = expression(theta), ylab = expression(kappa), main = "") # ensure axis labels are accurate points(a$theta1, a$theta2, lwd = 0.65) axis(side = 1, at = c(a$theta1hat, range(a$theta1))) axis(side = 2, at = c(a$theta2hat, range(a$theta2))) for (i in 1:nrow(pL)) { # ID & plot jump-center repair points in each quad m1R <- abs(a$theta1 - pR[i, 1]) / diff(range(a$theta1)) # match theta1 values m2R <- abs(a$theta2 - pR[i, 2]) / diff(range(a$theta2)) # match theta2 values m1L <- abs(a$theta1 - pL[i, 1]) / diff(range(a$theta1)) # match theta1 values m2L <- abs(a$theta2 - pL[i, 2]) / diff(range(a$theta2)) # match theta2 values fm <- which.min(m1R + m2R) + 1 # "right" points (nearest right boundary point) to <- which.min(m1L + m2L) - 1 # "to" points (nearest left boundary point) points(a$theta1[fm:to], a$theta2[fm:to], pch = 16, cex = 0.4, col = "blue") } points(pL, pch = 24, bg = "green") # points on left-side of inaccessible region points(pR, pch = 24, bg = "yellow") # points on right-side of inaccessible region points(pJ, pch = 24, bg = "red") # jump-center legend("topright", legend = c("CR boundary points", "jump-center (JC)", "JC left boundary", "JC right boundary", "JC repair points"), pch = c(1, rep(24, 3), 21), pt.bg = c("black", "red", "green", "yellow", "blue"), pt.cex = c(1, 1, 1, 1, 0.7), bty = "n", cex = 0.8) ## ----fig.width = 3.5, fig.height = 3.5, fig.show = 'hold', warning = FALSE, eval = FALSE---- # crplot(bc, alpha = 0.1, distn = "gamma", cen = cen, jumpshift = 0.1, showjump = TRUE) # #> [1] "Confidence region plot complete; made using 292 boundary points." ## ----fig.width = 3.5, fig.height = 3.5, fig.show = 'hold', echo = FALSE------- par(mar = c(4, 4.5, 2, 1.5)) # results of above chunk are input directly here to save processing time on vignette (otherwise prohibitively long for CRAN) a$theta1 <- c(11092.63, 10433.13, 9807.861, 8634.626, 8073.472, 7796.363, 7520.319, 6952.024, 6371.019, 5763.591, 5061.823, 4510.083, 4126.239, 4077.722, 4069.232, 4065.99, 4065.685, 4065.64, 4065.629, 4065.625, 3825.92, 3801.097, 3699.548, 3478.479, 3326.884, 3208.242, 3109.742, 3025.123, 2950.778, 2884.408, 2824.445, 2769.764, 2719.526, 2673.086, 2629.938, 2589.675, 2551.963, 2516.527, 2483.136, 2451.592, 2017.471, 1822.294, 1632.775, 1626.699, 1620.624, 1614.435, 1608.245, 1601.938, 1595.631, 1589.204, 1582.776, 1576.225, 1569.673, 1562.996, 1556.317, 1549.509, 1542.7, 1535.759, 1528.818, 1521.742, 1514.667, 1507.455, 1500.243, 1492.895, 1485.548, 1478.063, 1470.582, 1462.965, 1455.354, 1447.61, 1439.877, 1432.016, 1424.172, 1416.21, 1408.274, 1400.234, 1392.233, 1384.148, 1376.122, 1368.041, 1360.043, 1352.033, 1344.139, 1336.29, 1328.603, 1321.035, 1313.691, 1306.559, 1299.726, 1296.428, 1294.813, 1293.222, 1292.435, 1291.655, 1290.881, 1290.113, 1289.352, 1288.597, 1287.849, 1287.108, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.423, 1557.423, 1557.424, 1557.425, 1557.426, 1557.429, 1557.433, 1557.437, 1557.444, 1557.45, 1557.469, 1557.489, 1557.527, 1557.603, 1557.723, 1557.844, 1558.102, 1558.36, 1558.877, 1559.69, 1560.503, 1561.782, 1563.064, 1565.083, 1567.108, 1578.992, 1591.141, 1616.27, 1653.549, 1693.321, 1753.816, 1820.462, 2097.608, 2518.007, 2826.445, 3253.11, 3937.509, 4439.756, 5098.404, 5620.529, 6212.872, 6854.903, 7520.319, 8305.757, 9067.957, 9792.367, 10483.48, 11730.83, 12880.77, 14784.04, 16510.45, 19141.57, 21541.62, 25392.41, 27192.46, 28930.35, 35822.83, 42440.97, 48878.38, 49595.8, 50309.54, 51411.59, 51837.43, 52262.06, 52591.09, 52919.42, 53173.97, 53428.11, 53625.22, 53822.09, 53974.84, 54127.44, 54245.87, 54364.21, 54456.07, 54547.88, 54619.15, 54690.39, 54745.7, 54801, 54843.93, 54865.39, 54886.85, 54894.62, 54902.38, 54908.34, 54914.29, 54918.86, 54923.43, 54926.93, 54930.43, 54933.12, 54935.81, 54937.87, 54939.93, 54941.51, 54943.09, 54944.31, 54945.52, 373623.2, 376919.9, 379647, 380717.3, 381521.1, 381673.6, 381803.4, 381909, 381988.7, 382040.6, 382062.5, 382051.9, 382006.1, 381968.9, 381921.7, 381863.9, 381795.1, 381622.2, 381398.2, 380772.8, 379859.9, 370603.1, 290097.2, 255435.6, 237011.9, 223807.4, 213359.4, 204658.4, 197179.9, 190612, 184752.6, 179462.3, 174640.4, 170211.7, 166118.1, 162314, 158762.8, 155434.5, 152304.3, 149351.3, 146558.1, 143909.4, 141392.4, 138995.6, 136709.3, 107771.9, 84471.51, 79316.71, 74348.64, 64890.8, 60383.3, 55975.8, 51648.11, 47363.57, 46913.61, 46463.58, 45543.32, 44601.11, 43634.99, 42642.65, 41621.32, 40567.66, 39477.6, 38346.01, 37166.35, 35930.01, 34625.3, 33235.58, 31735.73, 30084.53, 28204.83, 25913.8, 23927.59, 23700.27, 23646.83, 23634.52, 23630.15, 23628.6, 21951.15, 21922.4, 21223.25, 18865.36, 17491.4, 16398.69, 15484.77, 13878.77, 12551.42) a$theta2 <- c(2.142689, 2.180192, 2.218991, 2.302236, 2.347989, 2.372294, 2.397765, 2.454672, 2.520183, 2.598519, 2.705241, 2.805313, 2.885938, 2.896894, 2.898831, 2.899572, 2.899642, 2.899652, 2.899655, 2.899655, 2.956911, 2.963133, 2.989204, 3.049657, 3.094367, 3.131423, 3.163691, 3.192584, 3.218923, 3.243234, 3.26588, 3.287123, 3.307159, 3.326142, 3.344191, 3.361404, 3.377863, 3.393633, 3.408773, 3.423331, 3.652807, 3.777007, 3.911765, 3.9163, 3.920846, 3.92549, 3.930144, 3.934898, 3.939663, 3.94453, 3.949408, 3.954388, 3.959379, 3.964476, 3.969581, 3.974793, 3.980012, 3.985339, 3.99067, 3.996109, 4.00155, 4.007096, 4.012641, 4.018288, 4.023928, 4.029665, 4.035387, 4.041197, 4.046983, 4.052846, 4.058671, 4.064555, 4.070382, 4.076246, 4.082028, 4.087813, 4.093484, 4.099114, 4.104586, 4.109956, 4.115112, 4.120087, 4.124772, 4.129175, 4.133189, 4.136795, 4.13989, 4.142422, 4.144296, 4.144962, 4.145222, 4.145432, 4.145517, 4.145589, 4.145648, 4.145694, 4.145726, 4.145744, 4.145748, 4.145737, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502057, 3.502056, 3.502055, 3.502054, 3.502052, 3.502049, 3.502043, 3.502038, 3.502029, 3.502021, 3.501996, 3.501971, 3.501922, 3.501822, 3.501666, 3.501509, 3.501173, 3.500838, 3.500166, 3.499112, 3.498058, 3.496403, 3.494747, 3.492148, 3.489547, 3.474448, 3.459287, 3.428771, 3.385445, 3.341559, 3.278935, 3.215022, 2.993232, 2.747052, 2.609885, 2.459118, 2.278387, 2.176874, 2.069931, 2.000331, 1.933303, 1.871587, 1.816841, 1.761516, 1.715315, 1.676807, 1.644061, 1.592815, 1.552607, 1.497032, 1.455406, 1.40337, 1.364588, 1.314333, 1.294595, 1.277305, 1.221507, 1.181113, 1.149826, 1.146713, 1.14368, 1.139119, 1.137395, 1.135696, 1.134393, 1.133105, 1.132114, 1.131132, 1.130375, 1.129623, 1.129042, 1.128464, 1.128018, 1.127573, 1.127228, 1.126885, 1.126619, 1.126353, 1.126147, 1.125942, 1.125783, 1.125703, 1.125624, 1.125595, 1.125566, 1.125544, 1.125522, 1.125505, 1.125488, 1.125475, 1.125462, 1.125453, 1.125443, 1.125435, 1.125427, 1.125421, 1.125416, 1.125411, 1.125407, 0.8893477, 0.8907084, 0.8925767, 0.8937651, 0.895167, 0.8955592, 0.8959693, 0.8963985, 0.8968482, 0.8973199, 0.8978152, 0.898336, 0.8988844, 0.8991703, 0.8994642, 0.8997662, 0.9000768, 0.9007257, 0.9014143, 0.9029291, 0.9046664, 0.9158979, 0.9809776, 1.010748, 1.028131, 1.041466, 1.05263, 1.062399, 1.071177, 1.079205, 1.086643, 1.093598, 1.100151, 1.106362, 1.112276, 1.117931, 1.123356, 1.128575, 1.133609, 1.138475, 1.143187, 1.147758, 1.152199, 1.156519, 1.160727, 1.223024, 1.291135, 1.309538, 1.32881, 1.370628, 1.393493, 1.418146, 1.444998, 1.474713, 1.478042, 1.481414, 1.488449, 1.495851, 1.503663, 1.511933, 1.520719, 1.530092, 1.540142, 1.550979, 1.562746, 1.575634, 1.589903, 1.605928, 1.624279, 1.6459, 1.672576, 1.708586, 1.743539, 1.747794, 1.748803, 1.749036, 1.749118, 1.749148, 1.782534, 1.783137, 1.79814, 1.854264, 1.891683, 1.924498, 1.954384, 2.013396, 2.069931) a$theta1hat <- 7520.318946 # MLE a$theta2hat <- 2.069931 # MLE pJ <- matrix(c(1632.775, 3.423331, 47363.57, 1.160727), nrow = 2, byrow = TRUE) # jump-center location pL <- matrix(c(1557.422, 3.502058, 21922.4, 1.783137), nrow = 2, byrow = TRUE) # jump-center left inaccessible region boundary pR <- matrix(c(4065.629, 2.899655, 54943.09, 1.125416), nrow = 2, byrow = TRUE) # jump-center right inaccessible region boundary #---> this point forward is universal for above inputs plot(c(a$theta1, a$theta1[1]), c(a$theta2, a$theta2[1]), type = "l", axes = FALSE, xlab = expression(theta), ylab = expression(kappa), main = "") # ensure axis labels are accurate points(a$theta1, a$theta2, lwd = 0.65) axis(side = 1, at = c(a$theta1hat, range(a$theta1))) axis(side = 2, at = c(a$theta2hat, range(a$theta2))) for (i in 1:nrow(pL)) { # ID & plot jump-center repair points in each quad m1R <- abs(a$theta1 - pR[i, 1]) / diff(range(a$theta1)) # match theta1 values m2R <- abs(a$theta2 - pR[i, 2]) / diff(range(a$theta2)) # match theta2 values m1L <- abs(a$theta1 - pL[i, 1]) / diff(range(a$theta1)) # match theta1 values m2L <- abs(a$theta2 - pL[i, 2]) / diff(range(a$theta2)) # match theta2 values fm <- which.min(m1R + m2R) + 1 # "right" points (nearest right boundary point) to <- which.min(m1L + m2L) - 1 # "to" points (nearest left boundary point) points(a$theta1[fm:to], a$theta2[fm:to], pch = 16, cex = 0.4, col = "blue") } points(pL, pch = 24, bg = "green") # points on left-side of inaccessible region points(pR, pch = 24, bg = "yellow") # points on right-side of inaccessible region points(pJ, pch = 24, bg = "red") # jump-center legend("topright", legend = c("CR boundary points", "jump-center (JC)", "JC left boundary", "JC right boundary", "JC repair points"), pch = c(1, rep(24, 3), 21), pt.bg = c("black", "red", "green", "yellow", "blue"), pt.cex = c(1, 1, 1, 1, 0.7), bty = "n", cex = 0.8) ## ----fig.width = 3.5, fig.height = 3.5, fig.show = 'hold', warning = FALSE, eval = FALSE---- # crplot(bc, alpha = 0.1, distn = "gamma", cen = cen, showjump = TRUE, ellipse_n = 80) # #> [1] "Confidence region plot complete; made using 387 boundary points." ## ----fig.width = 3.5, fig.height = 3.5, fig.show = 'hold', echo = FALSE------- par(mar = c(4, 4.5, 2, 1.5)) # results of above chunk are input directly here to save processing time on vignette (otherwise prohibitively long for CRAN) a$theta1 <- c(12444.4, 11782.27, 11241.81, 10828.5, 10494.14, 10212.32, 9967.056, 9748.017, 9548.122, 9362.297, 9186.756, 9018.561, 8855.346, 8695.125, 8536.167, 8376.885, 8215.761, 8051.265, 7881.779, 7705.499, 7520.319, 7312.922, 7062.707, 6748.991, 6330.569, 5698.124, 4793.889, 4163.64, 4105.016, 4078.244, 4067.186, 4066.112, 4065.716, 4065.658, 4065.636, 4065.628, 4065.625, 4065.625, 4065.625, 3854.732, 3636.822, 3498.94, 3400.909, 3326.124, 3266.188, 3216.304, 3173.494, 3135.783, 3101.797, 3067.036, 3041.216, 3013.248, 2986.119, 2959.39, 2932.655, 2905.526, 2877.607, 2848.473, 2817.651, 2784.585, 2748.601, 2708.847, 2664.21, 2613.173, 2553.59, 2482.254, 2394.059, 2280.098, 2122.504, 1875.169, 1585.041, 1432.731, 1359.298, 1326.368, 1311.584, 1298.246, 1298.138, 1298.031, 1297.815, 1297.387, 1296.536, 1294.858, 1293.213, 1292.403, 1291.602, 1291.204, 1290.809, 1290.416, 1290.025, 1289.83, 1289.636, 1289.442, 1289.249, 1289.057, 1288.865, 1288.673, 1288.482, 1288.292, 1288.102, 1287.913, 1287.725, 1287.536, 1287.349, 1287.162, 1286.975, 1286.79, 1286.604, 1286.42, 1286.235, 1286.052, 1285.869, 1285.504, 1285.142, 1284.782, 1284.069, 1283.365, 1282.671, 1281.31, 1279.986, 1277.453, 1275.076, 1261.929, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.423, 1557.423, 1557.424, 1557.427, 1557.429, 1557.433, 1557.443, 1557.457, 1557.472, 1557.495, 1557.518, 1557.61, 1557.701, 1557.885, 1558.253, 1558.989, 1560.142, 1561.297, 1563.107, 1564.923, 1567.774, 1570.64, 1575.15, 1579.697, 1711.66, 1878.653, 2192.796, 2512.726, 2830.903, 3140.541, 3435.873, 3712.651, 3968.442, 4202.573, 4415.78, 4609.738, 4786.617, 4948.747, 5098.404, 5237.701, 5368.533, 5492.576, 5611.299, 5725.99, 5837.779, 5947.666, 6056.54, 6165.204, 6274.384, 6384.749, 6496.914, 6611.452, 6728.898, 6849.753, 6974.483, 7103.522, 7237.272, 7376.095, 7520.319, 7678.279, 7861.514, 8076.275, 8330.79, 8635.963, 9006.358, 9461.527, 10027.88, 10741.38, 11651.62, 12828.66, 14375.38, 16451.21, 17927.17, 19321.31, 23466.6, 29862.82, 35540.48, 40828.61, 42604.55, 44351.76, 47287.44, 49747.96, 51818.67, 52695.07, 53566.51, 53736.02, 53905.34, 54035.88, 54166.31, 54266.89, 54367.41, 54444.93, 54522.42, 54582.19, 54641.93, 54688.02, 54734.09, 54769.64, 54805.18, 54832.6, 54860.01, 54881.16, 54902.31, 54918.62, 54926.78, 54934.94, 54937.87, 54940.79, 54943.04, 54944.16, 54944.73, 54945.29, 54945.38, 54945.47, 54945.55, 54945.58, 54945.62, 54945.63, 54945.64, 54945.65, 54945.66, 54945.67, 54945.68, 54945.68, 55577.95, 60732.38, 66451.55, 72747.55, 79640.83, 87157.84, 95329.62, 104190.8, 113779.1, 124134.7, 135300, 147319.4, 160238.4, 174103.6, 188961.1, 204856.1, 221829.6, 239915.6, 259134.3, 279480.1, 300897.6, 323226.4, 334878.7, 340529.5, 346050, 355142.6, 359496.1, 363693.9, 364261.5, 364543.9, 364825.4, 364949.3, 365073.1, 365157, 365198.9, 365219.9, 365240.9, 365246.8, 365252.7, 365256.7, 365258.7, 365259.7, 365260.7, 365261, 365261.3, 365261.5, 365261.6, 365261.7, 365261.7, 365261.7, 365261.7, 365261.7, 365261.7, 365261.7, 365261.7, 191047.1, 191047.1, 191047, 191046.9, 191046.4, 191044.2, 191029.6, 190984.5, 190682.1, 189774.3, 184504.2, 168701.4, 141255.2, 116065.6, 100938.6, 89979.57, 81350.09, 66143.49, 59459.21, 53242.26, 49380.11, 45660.38, 42135.91, 40032.02, 38605.07, 37556.79, 36742.69, 36083.81, 35533.19, 35061.06, 34647.63, 34279.21, 33945.99, 33640.79, 33358.21, 33094.1, 32845.23, 32609.04, 32383.49, 32166.89, 31947.07, 31711.08, 31455.86, 31177.46, 30870.76, 30528.94, 30142.71, 29699.09, 29179.3, 28555, 27781.09, 26780.32, 25405.32, 23318.49, 21959.33, 21959.28, 21959.1, 21957.35, 21944.69, 21892.63, 21583.42, 18649.03, 17079.06, 16094.08, 15251.3, 13825.24, 13106.33) a$theta2 <- c(2.074855, 2.106707, 2.134641, 2.15732, 2.176582, 2.193508, 2.208785, 2.222888, 2.236157, 2.248851, 2.261176, 2.273304, 2.285385, 2.297556, 2.309952, 2.322707, 2.335967, 2.349893, 2.364673, 2.380534, 2.397765, 2.417798, 2.443077, 2.476644, 2.525061, 2.607669, 2.751771, 2.877616, 2.890708, 2.896776, 2.899299, 2.899544, 2.899635, 2.899648, 2.899653, 2.899655, 2.899655, 2.899655, 2.899655, 2.949761, 3.005823, 3.043835, 3.072181, 3.094598, 3.113086, 3.128843, 3.142645, 3.155021, 3.166357, 3.178133, 3.187003, 3.19673, 3.206288, 3.215825, 3.225485, 3.235415, 3.24577, 3.256726, 3.268489, 3.281308, 3.295502, 3.311484, 3.329822, 3.351314, 3.377146, 3.409177, 3.450546, 3.507128, 3.591797, 3.741906, 3.947688, 4.064021, 4.115583, 4.134293, 4.140692, 4.144616, 4.144638, 4.14466, 4.144703, 4.144786, 4.144943, 4.145215, 4.145433, 4.14552, 4.145594, 4.145625, 4.145653, 4.145678, 4.145698, 4.145707, 4.145715, 4.145723, 4.145729, 4.145734, 4.145739, 4.145743, 4.145745, 4.145747, 4.145748, 4.145748, 4.145747, 4.145745, 4.145742, 4.145739, 4.145734, 4.145728, 4.145722, 4.145715, 4.145706, 4.145697, 4.145687, 4.145664, 4.145637, 4.145606, 4.145534, 4.145447, 4.145344, 4.145093, 4.144781, 4.143972, 4.14291, 4.125071, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502057, 3.502056, 3.502055, 3.502052, 3.502049, 3.502043, 3.502031, 3.502012, 3.501993, 3.501963, 3.501933, 3.501814, 3.501694, 3.501455, 3.500977, 3.500021, 3.498526, 3.497031, 3.494693, 3.492353, 3.488694, 3.485031, 3.4793, 3.47356, 3.322073, 3.163035, 2.929715, 2.749667, 2.608102, 2.495359, 2.404747, 2.331371, 2.271503, 2.222233, 2.18127, 2.146817, 2.117465, 2.092121, 2.069931, 2.050235, 2.032515, 2.016366, 2.001467, 1.987562, 1.974445, 1.96195, 1.94994, 1.938302, 1.926944, 1.915788, 1.904769, 1.893833, 1.882937, 1.872044, 1.861125, 1.850159, 1.839129, 1.828024, 1.816841, 1.804984, 1.791714, 1.776781, 1.759891, 1.740699, 1.718815, 1.693811, 1.665247, 1.632701, 1.595811, 1.554307, 1.507998, 1.456722, 1.425937, 1.400215, 1.337915, 1.268653, 1.223472, 1.190059, 1.180233, 1.171135, 1.156974, 1.146061, 1.13747, 1.133984, 1.1306, 1.129951, 1.129306, 1.128811, 1.128318, 1.127938, 1.127561, 1.12727, 1.12698, 1.126756, 1.126534, 1.126362, 1.126191, 1.126058, 1.125926, 1.125825, 1.125723, 1.125645, 1.125567, 1.125506, 1.125476, 1.125446, 1.125435, 1.125424, 1.125416, 1.125412, 1.12541, 1.125408, 1.125407, 1.125407, 1.125407, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.123088, 1.105524, 1.088414, 1.071897, 1.056046, 1.040893, 1.026447, 1.012702, 0.9996437, 0.9872559, 0.9755203, 0.9644199, 0.9539401, 0.9440703, 0.9348052, 0.9261468, 0.9181069, 0.9107118, 0.9040086, 0.8980787, 0.8930627, 0.8892178, 0.8878554, 0.8873871, 0.8870733, 0.8869384, 0.8870906, 0.8874169, 0.8874775, 0.8875092, 0.887542, 0.8875568, 0.8875718, 0.8875821, 0.8875873, 0.8875899, 0.8875925, 0.8875933, 0.887594, 0.8875945, 0.8875947, 0.8875949, 0.887595, 0.887595, 0.8875951, 0.8875951, 0.8875951, 0.8875951, 0.8875951, 0.8875951, 0.8875951, 0.8875951, 0.8875951, 0.8875951, 0.8875951, 1.078664, 1.078664, 1.078664, 1.078664, 1.078665, 1.078667, 1.078685, 1.078742, 1.079118, 1.080252, 1.086964, 1.108524, 1.152444, 1.203193, 1.240881, 1.273014, 1.302098, 1.364643, 1.39846, 1.434769, 1.460298, 1.487543, 1.516256, 1.534984, 1.54846, 1.558796, 1.567096, 1.573998, 1.579898, 1.585057, 1.589653, 1.593811, 1.597625, 1.601163, 1.604478, 1.60761, 1.610594, 1.613454, 1.616211, 1.618884, 1.621622, 1.62459, 1.627834, 1.631414, 1.635408, 1.639924, 1.64511, 1.651181, 1.658453, 1.667425, 1.678928, 1.69447, 1.717174, 1.755069, 1.782363, 1.782364, 1.782368, 1.782405, 1.78267, 1.783763, 1.790328, 1.859899, 1.903721, 1.934195, 1.962417, 2.015527, 2.045317) a$theta1hat <- 7520.318946 # MLE a$theta2hat <- 2.069931 # MLE pJ <- matrix(c(1875.169, 3.176941, 32166.89, 1.252914), nrow = 2, byrow = TRUE) # jump-center location pL <- matrix(c(1557.422, 3.502058, 21959.28, 1.782364), nrow = 2, byrow = TRUE) # jump-center left inaccessible region boundary pR <- matrix(c(4065.625, 2.899655, 54945.67, 1.125406), nrow = 2, byrow = TRUE) # jump-center right inaccessible region boundary #---> this point forward is universal for above inputs plot(c(a$theta1, a$theta1[1]), c(a$theta2, a$theta2[1]), type = "l", axes = FALSE, xlab = expression(theta), ylab = expression(kappa), main = "") # ensure axis labels are accurate points(a$theta1, a$theta2, lwd = 0.65) axis(side = 1, at = c(a$theta1hat, range(a$theta1))) axis(side = 2, at = c(a$theta2hat, range(a$theta2))) for (i in 1:nrow(pL)) { # ID & plot jump-center repair points in each quad m1R <- abs(a$theta1 - pR[i, 1]) / diff(range(a$theta1)) # match theta1 values m2R <- abs(a$theta2 - pR[i, 2]) / diff(range(a$theta2)) # match theta2 values m1L <- abs(a$theta1 - pL[i, 1]) / diff(range(a$theta1)) # match theta1 values m2L <- abs(a$theta2 - pL[i, 2]) / diff(range(a$theta2)) # match theta2 values fm <- which.min(m1R + m2R) + 1 # "right" points (nearest right boundary point) to <- which.min(m1L + m2L) - 1 # "to" points (nearest left boundary point) points(a$theta1[fm:to], a$theta2[fm:to], pch = 16, cex = 0.4, col = "blue") } points(pL, pch = 24, bg = "green") # points on left-side of inaccessible region points(pR, pch = 24, bg = "yellow") # points on right-side of inaccessible region points(pJ, pch = 24, bg = "red") # jump-center legend("topright", legend = c("CR boundary points", "jump-center (JC)", "JC left boundary", "JC right boundary", "JC repair points"), pch = c(1, rep(24, 3), 21), pt.bg = c("black", "red", "green", "yellow", "blue"), pt.cex = c(1, 1, 1, 1, 0.7), bty = "n", cex = 0.8) ## ----fig.width = 3.5, fig.height = 3.5, fig.show = 'hold', warning = FALSE, eval = FALSE---- # crplot(bc, alpha = 0.1, distn = "gamma", cen = cen, maxcount = 50) # #> [1] "Confidence region plot complete; made using 476 boundary points." ## ----fig.width = 3.5, fig.height = 3.5, fig.show = 'hold', echo = FALSE------- par(mar = c(4, 4.5, 2, 1.5)) # results of above chunk are input directly here to save processing time on vignette (otherwise prohibitively long for CRAN) a$theta1 <- c(11092.63, 10433.13, 9807.861, 8634.626, 8073.472, 7796.363, 7520.319, 6952.024, 6371.019, 5763.591, 5061.823, 4510.083, 4126.239, 4077.722, 4069.232, 4065.99, 4065.685, 4065.64, 4065.629, 4065.625, 4065.625, 4065.625, 4065.625, 3988.343, 3894.444, 3810.222, 3733.785, 3663.773, 3599.172, 3539.201, 3483.245, 3430.811, 3381.495, 3334.96, 3290.926, 3249.154, 3209.438, 3171.602, 3135.491, 3100.97, 3067.919, 2422.231, 2142.175, 1874.155, 1867.84, 1861.525, 1855.117, 1848.708, 1842.202, 1835.696, 1829.091, 1822.484, 1815.776, 1809.065, 1802.25, 1795.433, 1788.508, 1781.579, 1774.54, 1767.498, 1760.341, 1753.179, 1745.9, 1738.616, 1731.211, 1723.8, 1716.263, 1708.72, 1701.048, 1693.368, 1685.555, 1677.734, 1669.775, 1661.808, 1653.698, 1645.579, 1637.313, 1629.038, 1620.61, 1612.173, 1603.58, 1594.977, 1586.214, 1577.441, 1568.504, 1559.558, 1550.445, 1541.325, 1532.036, 1522.743, 1513.281, 1503.819, 1494.192, 1484.571, 1474.791, 1465.028, 1455.117, 1445.238, 1435.231, 1425.277, 1415.224, 1405.253, 1395.227, 1385.323, 1375.423, 1365.698, 1356.056, 1346.657, 1337.441, 1328.549, 1319.957, 1311.779, 1304.025, 1300.335, 1296.776, 1295.907, 1295.047, 1293.352, 1292.518, 1291.693, 1291.284, 1290.877, 1290.473, 1290.07, 1289.87, 1289.67, 1289.471, 1289.273, 1289.075, 1288.877, 1288.68, 1288.484, 1288.289, 1288.093, 1287.996, 1287.899, 1287.802, 1287.705, 1287.512, 1287.319, 1287.127, 1286.935, 1286.744, 1286.554, 1286.364, 1286.175, 1285.986, 1285.798, 1285.424, 1285.052, 1284.683, 1284.316, 1284.133, 1283.951, 1283.232, 1282.522, 1281.822, 1281.132, 1278.471, 1275.971, 1274.783, 1273.636, 1269.503, 1266.048, 1261.273, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.423, 1557.423, 1557.424, 1557.425, 1557.426, 1557.429, 1557.433, 1557.437, 1557.444, 1557.45, 1557.469, 1557.489, 1557.527, 1557.603, 1557.723, 1557.844, 1558.102, 1558.36, 1558.877, 1559.69, 1560.503, 1561.782, 1563.064, 1565.083, 1567.108, 1578.992, 1591.141, 1616.27, 1653.549, 1693.321, 1753.816, 1820.462, 2097.608, 2518.007, 2826.445, 3253.11, 3937.509, 4439.756, 5098.404, 5620.529, 6212.872, 6854.903, 7520.319, 8305.757, 9067.957, 9792.367, 10483.48, 11730.83, 12880.77, 14784.04, 16510.45, 19141.57, 21541.62, 25392.41, 27192.46, 28930.35, 35822.83, 42440.97, 48878.38, 49595.8, 50309.54, 51411.59, 51837.43, 52262.06, 52591.09, 52919.42, 53173.97, 53428.11, 53625.22, 53822.09, 53974.84, 54127.44, 54245.87, 54364.21, 54456.07, 54547.88, 54619.15, 54690.39, 54745.7, 54801, 54843.93, 54865.39, 54886.85, 54894.62, 54902.38, 54908.34, 54914.29, 54918.86, 54923.43, 54926.93, 54930.43, 54933.12, 54935.81, 54937.87, 54939.93, 54941.51, 54943.09, 54944.31, 54944.91, 54945.22, 54945.52, 54945.57, 54945.59, 54945.62, 54945.63, 54945.64, 54945.65, 54945.65, 54945.66, 54945.66, 54945.67, 54945.67, 54945.67, 54945.68, 54945.68, 54945.68, 54945.68, 54945.68, 54945.68, 54945.68, 54945.68, 54945.68, 54945.68, 54945.68, 54945.68, 54945.68, 54945.68, 54945.68, 54945.68, 54945.68, 54945.68, 54945.68, 54945.68, 54945.68, 271649.9, 274925.9, 278222, 281537.7, 284872.6, 288226.3, 291598.2, 294987.8, 298394.3, 301816.9, 305254.8, 308706.8, 312172, 315648.9, 319136.1, 322631.8, 326134.1, 329640.5, 333148.6, 336655.1, 340156.5, 343648.3, 347125.6, 350582, 354010, 355715, 356560.2, 356980.9, 357400.3, 357776.8, 358152.3, 358408.8, 358664.7, 358839.7, 359014.4, 359133.9, 359193.6, 359253.2, 359279.1, 359305, 359322.7, 359331.6, 359340.4, 359344.2, 359346.2, 359347.1, 359347.6, 359347.8, 359348.1, 359348.1, 359348.1, 359348.2, 359348.2, 359348.2, 359348.2, 359348.2, 359348.2, 359348.2, 359348.2, 359348.2, 359348.2, 359348.2, 183812.8, 183812.8, 183812.7, 183812.6, 183811.3, 183807.3, 183794.6, 183754.6, 183571.5, 182378.2, 169902.3, 158854.7, 150907.6, 144566.4, 139242.6, 134633.4, 130559.2, 126903.5, 123585.6, 120547, 117744.1, 115142.7, 112716.3, 110443.2, 108305.6, 106289, 104380.9, 102570.8, 100849.8, 99209.97, 97644.63, 96147.8, 94714.22, 93339.23, 82395.09, 72829.27, 68445.41, 64288.49, 56554.39, 52979.94, 49542.88, 43020.31, 39922.33, 36897.95, 33932.67, 30999.44, 30764, 30528.54, 30049.42, 29561.55, 29064.33, 28557.06, 28038.98, 27509.18, 26966.64, 26410.11, 25838.18, 25249.09, 24640.76, 24010.57, 23355.24, 22670.53, 21959.33, 21959.33, 21959.33, 21959.25, 21958.02, 21951.15, 21922.4, 21223.25, 18865.36, 17491.4, 16398.69, 15484.77, 13878.77, 12551.42) a$theta2 <- c(2.142689, 2.180192, 2.218991, 2.302236, 2.347989, 2.372294, 2.397765, 2.454672, 2.520183, 2.598519, 2.705241, 2.805313, 2.885938, 2.896894, 2.898831, 2.899572, 2.899642, 2.899652, 2.899655, 2.899655, 2.899655, 2.899655, 2.899655, 2.917576, 2.94003, 2.960839, 2.980301, 2.998633, 3.015997, 3.032518, 3.048297, 3.063413, 3.077933, 3.091912, 3.105397, 3.118428, 3.131039, 3.14326, 3.155118, 3.166635, 3.177832, 3.437111, 3.580792, 3.742569, 3.746706, 3.750858, 3.755088, 3.759335, 3.763662, 3.768006, 3.772433, 3.776879, 3.781411, 3.785962, 3.790603, 3.795263, 3.800016, 3.80479, 3.80966, 3.814552, 3.819543, 3.824558, 3.829675, 3.834817, 3.840066, 3.845339, 3.850724, 3.856135, 3.86166, 3.867213, 3.872885, 3.878584, 3.884407, 3.890258, 3.896237, 3.902245, 3.908384, 3.914553, 3.920856, 3.927189, 3.93366, 3.940158, 3.946798, 3.953464, 3.960272, 3.967103, 3.974076, 3.981067, 3.988198, 3.99534, 4.002616, 4.009892, 4.017292, 4.024677, 4.032169, 4.039626, 4.047163, 4.054636, 4.062153, 4.069565, 4.076967, 4.084211, 4.091373, 4.098303, 4.105056, 4.111485, 4.117614, 4.123303, 4.128546, 4.133216, 4.137276, 4.140619, 4.143188, 4.144155, 4.1449, 4.145051, 4.145187, 4.145416, 4.145509, 4.145586, 4.145619, 4.145649, 4.145674, 4.145696, 4.145706, 4.145714, 4.145722, 4.145728, 4.145734, 4.145739, 4.145742, 4.145745, 4.145747, 4.145748, 4.145748, 4.145748, 4.145747, 4.145747, 4.145745, 4.145742, 4.145738, 4.145733, 4.145727, 4.14572, 4.145712, 4.145703, 4.145693, 4.145683, 4.145658, 4.14563, 4.145597, 4.145561, 4.145541, 4.14552, 4.145428, 4.14532, 4.145196, 4.145056, 4.144333, 4.143348, 4.142756, 4.142097, 4.138816, 4.134455, 4.122628, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502057, 3.502056, 3.502055, 3.502054, 3.502052, 3.502049, 3.502043, 3.502038, 3.502029, 3.502021, 3.501996, 3.501971, 3.501922, 3.501822, 3.501666, 3.501509, 3.501173, 3.500838, 3.500166, 3.499112, 3.498058, 3.496403, 3.494747, 3.492148, 3.489547, 3.474448, 3.459287, 3.428771, 3.385445, 3.341559, 3.278935, 3.215022, 2.993232, 2.747052, 2.609885, 2.459118, 2.278387, 2.176874, 2.069931, 2.000331, 1.933303, 1.871587, 1.816841, 1.761516, 1.715315, 1.676807, 1.644061, 1.592815, 1.552607, 1.497032, 1.455406, 1.40337, 1.364588, 1.314333, 1.294595, 1.277305, 1.221507, 1.181113, 1.149826, 1.146713, 1.14368, 1.139119, 1.137395, 1.135696, 1.134393, 1.133105, 1.132114, 1.131132, 1.130375, 1.129623, 1.129042, 1.128464, 1.128018, 1.127573, 1.127228, 1.126885, 1.126619, 1.126353, 1.126147, 1.125942, 1.125783, 1.125703, 1.125624, 1.125595, 1.125566, 1.125544, 1.125522, 1.125505, 1.125488, 1.125475, 1.125462, 1.125453, 1.125443, 1.125435, 1.125427, 1.125421, 1.125416, 1.125411, 1.125409, 1.125408, 1.125407, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 1.125406, 0.9002247, 0.8993064, 0.8984122, 0.8975426, 0.8966981, 0.8958793, 0.8950869, 0.8943215, 0.8935841, 0.8928756, 0.892197, 0.8915496, 0.8909346, 0.8903535, 0.889808, 0.8893002, 0.8888321, 0.8884062, 0.8880256, 0.8876934, 0.8874137, 0.8871909, 0.8870306, 0.8869393, 0.8869247, 0.8869492, 0.8869699, 0.8869825, 0.8869965, 0.8870104, 0.8870256, 0.8870367, 0.8870483, 0.8870567, 0.8870653, 0.8870714, 0.8870744, 0.8870776, 0.8870789, 0.8870803, 0.8870812, 0.8870817, 0.8870822, 0.8870824, 0.8870825, 0.8870825, 0.8870826, 0.8870826, 0.8870826, 0.8870826, 0.8870826, 0.8870826, 0.8870826, 0.8870826, 0.8870826, 0.8870826, 0.8870826, 0.8870826, 0.8870826, 0.8870826, 0.8870826, 0.8870826, 1.087861, 1.087861, 1.087861, 1.087861, 1.087863, 1.087868, 1.087885, 1.087937, 1.088175, 1.089735, 1.106803, 1.123214, 1.135896, 1.146615, 1.15607, 1.164624, 1.172493, 1.179819, 1.1867, 1.193208, 1.199397, 1.205309, 1.210978, 1.216431, 1.22169, 1.226774, 1.231699, 1.236478, 1.241123, 1.245644, 1.250049, 1.254348, 1.258546, 1.26265, 1.298367, 1.335043, 1.354025, 1.37356, 1.414765, 1.436424, 1.459169, 1.508755, 1.535997, 1.565494, 1.597778, 1.633725, 1.636811, 1.639929, 1.646377, 1.653088, 1.660086, 1.667395, 1.675047, 1.683074, 1.691517, 1.700423, 1.709848, 1.719861, 1.730544, 1.742, 1.754361, 1.767799, 1.782363, 1.782363, 1.782363, 1.782365, 1.782391, 1.782534, 1.783137, 1.79814, 1.854264, 1.891683, 1.924498, 1.954384, 2.013396, 2.069931) a$theta1hat <- 7520.318946 # MLE a$theta2hat <- 2.069931 # MLE #---> this point forward is universal for above inputs plot(c(a$theta1, a$theta1[1]), c(a$theta2, a$theta2[1]), type = "l", axes = FALSE, xlab = expression(theta), ylab = expression(kappa), main = "") # ensure axis labels are accurate points(a$theta1, a$theta2, lwd = 0.65) axis(side = 1, at = c(a$theta1hat, range(a$theta1))) axis(side = 2, at = c(a$theta2hat, range(a$theta2))) ## ----fig.show = 'hold', warning = FALSE, eval = FALSE------------------------- # x <- crplot(bc, alpha = 0.1, distn = "gamma", cen = cen, showjump = TRUE, # jumpshift = 0.07, jumpuphill = 0.05, jumpinfo = TRUE) # #> [1] "Confidence region plot complete; made using 240 boundary points." # plot(c(x$theta, x$theta[1]), c(x$kappa, x$kappa[1]), type = "l", axes = FALSE, # xlab = expression(theta), ylab = expression(kappa)) # axis(side = 1, at = range(x$theta)) # axis(side = 2, at = range(x$kappa)) ## ----fig.show = 'hold', echo = FALSE------------------------------------------ par(mar = c(4, 4.5, 2, 1.5)) # results of above chunk are input directly here to save processing time on vignette (otherwise prohibitively long for CRAN) a$theta1 <- c(11092.63, 10433.13, 9807.861, 8634.626, 8073.472, 7796.363, 7520.319, 6952.024, 6371.019, 5763.591, 5061.823, 4510.083, 4126.239, 4077.722, 4069.232, 4065.99, 4065.685, 4065.64, 4065.629, 4065.625, 1483.341, 1499.991, 1517.076, 1534.788, 1552.943, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.423, 1557.423, 1557.424, 1557.425, 1557.426, 1557.429, 1557.433, 1557.437, 1557.444, 1557.45, 1557.469, 1557.489, 1557.527, 1557.603, 1557.723, 1557.844, 1558.102, 1558.36, 1558.877, 1559.69, 1560.503, 1561.782, 1563.064, 1565.083, 1567.108, 1578.992, 1591.141, 1616.27, 1653.549, 1693.321, 1753.816, 1820.462, 2097.608, 2518.007, 2826.445, 3253.11, 3937.509, 4439.756, 5098.404, 5620.529, 6212.872, 6854.903, 7520.319, 8305.757, 9067.957, 9792.367, 10483.48, 11730.83, 12880.77, 14784.04, 16510.45, 19141.57, 21541.62, 25392.41, 27192.46, 28930.35, 35822.83, 42440.97, 48878.38, 49595.8, 50309.54, 51411.59, 51837.43, 52262.06, 52591.09, 52919.42, 53173.97, 53428.11, 53625.22, 53822.09, 53974.84, 54127.44, 54245.87, 54364.21, 54456.07, 54547.88, 54619.15, 54690.39, 54745.7, 54801, 54843.93, 54865.39, 54886.85, 54894.62, 54902.38, 54908.34, 54914.29, 54918.86, 54923.43, 54926.93, 54930.43, 54933.12, 54935.81, 54937.87, 54939.93, 54941.51, 54943.09, 54944.31, 54945.52, 57217.34, 63771.85, 72528.15, 82785.21, 94037.64, 106050.1, 118733.2, 132062, 146038.9, 160677.3, 175994.4, 192006.8, 208727.9, 226164.7, 244313.4, 263153.7, 282636.9, 302665.4, 323049.4, 343404.4, 362852.3, 371710.1, 375523.5, 378740.9, 380484.8, 381161.3, 381440.1, 381674.1, 381772.9, 381858.6, 381896.4, 381930.6, 381961.1, 381987.9, 382010.8, 382020.8, 382029.7, 382037.6, 382044.5, 382050.3, 382055.1, 382058.8, 382061.3, 382062.7, 382063, 382062.1, 382060, 382056.8, 382052.3, 382046.5, 382039.5, 382031.2, 382021.6, 382010.7, 381998.4, 381969.6, 381935.1, 381894.6, 381795, 381668.6, 381513.2, 381326.3, 381220.1, 381163.6, 381104.8, 379695.7, 377165.9, 281718.9, 224935.4, 181174.9, 158666.2, 143303.7, 131713.7, 117615.1, 105104.6, 93802.35, 83459.78, 78636.93, 73979.37, 65086.54, 60833.62, 56668.29, 52572.15, 48513.22, 46299.82, 44078.33, 38943.62, 35777.59, 32403.98, 29268.35, 25363.77, 23505.32, 23144.91, 23109.1, 23100.73, 23098.79, 23098.34, 23098.17, 23098.14, 23098.12, 23098.12, 21951.15, 21922.4, 21223.25, 18865.36, 17491.4, 16398.69, 15484.77, 13878.77, 12551.42) a$theta2 <- c(2.142689, 2.180192, 2.218991, 2.302236, 2.347989, 2.372294, 2.397765, 2.454672, 2.520183, 2.598519, 2.705241, 2.805313, 2.885938, 2.896894, 2.898831, 2.899572, 2.899642, 2.899652, 2.899655, 2.899655, 3.604532, 3.580366, 3.556287, 3.532044, 3.507908, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502057, 3.502056, 3.502055, 3.502054, 3.502052, 3.502049, 3.502043, 3.502038, 3.502029, 3.502021, 3.501996, 3.501971, 3.501922, 3.501822, 3.501666, 3.501509, 3.501173, 3.500838, 3.500166, 3.499112, 3.498058, 3.496403, 3.494747, 3.492148, 3.489547, 3.474448, 3.459287, 3.428771, 3.385445, 3.341559, 3.278935, 3.215022, 2.993232, 2.747052, 2.609885, 2.459118, 2.278387, 2.176874, 2.069931, 2.000331, 1.933303, 1.871587, 1.816841, 1.761516, 1.715315, 1.676807, 1.644061, 1.592815, 1.552607, 1.497032, 1.455406, 1.40337, 1.364588, 1.314333, 1.294595, 1.277305, 1.221507, 1.181113, 1.149826, 1.146713, 1.14368, 1.139119, 1.137395, 1.135696, 1.134393, 1.133105, 1.132114, 1.131132, 1.130375, 1.129623, 1.129042, 1.128464, 1.128018, 1.127573, 1.127228, 1.126885, 1.126619, 1.126353, 1.126147, 1.125942, 1.125783, 1.125703, 1.125624, 1.125595, 1.125566, 1.125544, 1.125522, 1.125505, 1.125488, 1.125475, 1.125462, 1.125453, 1.125443, 1.125435, 1.125427, 1.125421, 1.125416, 1.125411, 1.125407, 1.117252, 1.096153, 1.072437, 1.049464, 1.028608, 1.010034, 0.993517, 0.9787697, 0.9655359, 0.9536073, 0.942821, 0.9330526, 0.9242103, 0.9162307, 0.909078, 0.9027452, 0.8972609, 0.892705, 0.8892422, 0.8872046, 0.8873349, 0.8887878, 0.8900543, 0.8918315, 0.8934649, 0.8944449, 0.8949842, 0.8955606, 0.8958643, 0.8961787, 0.8963402, 0.8965045, 0.8966719, 0.8968423, 0.897016, 0.897104, 0.8971929, 0.8972826, 0.8973731, 0.8974645, 0.8975568, 0.89765, 0.8977441, 0.8978391, 0.8979351, 0.898032, 0.8981299, 0.8982287, 0.8983286, 0.8984295, 0.8985314, 0.8986344, 0.8987384, 0.8988436, 0.8989498, 0.8991657, 0.8993862, 0.8996115, 0.9000775, 0.9005651, 0.9010762, 0.9016128, 0.9018913, 0.9020333, 0.9021772, 0.9049458, 0.9086121, 0.9878909, 1.040294, 1.09132, 1.123505, 1.148818, 1.170232, 1.199686, 1.229817, 1.261259, 1.29463, 1.31208, 1.330309, 1.369683, 1.391109, 1.414105, 1.439018, 1.466397, 1.482652, 1.50005, 1.545203, 1.577264, 1.61596, 1.657195, 1.717887, 1.751489, 1.758431, 1.759129, 1.759292, 1.75933, 1.759339, 1.759342, 1.759343, 1.759343, 1.759343, 1.782534, 1.783137, 1.79814, 1.854264, 1.891683, 1.924498, 1.954384, 2.013396, 2.069931) a$theta1hat <- 7520.318946 # MLE a$theta2hat <- 2.069931 # MLE pJ <- matrix(c(4627.62, 2.743556, 48513.22, 1.170232), nrow = 2, byrow = TRUE) # jump-center location pL <- matrix(c(1557.422, 3.502058, 21922.4, 1.783137), nrow = 2, byrow = TRUE) # jump-center left inaccessible region boundary pR <- matrix(c(4065.629, 2.899655, 54943.09, 1.125416), nrow = 2, byrow = TRUE) # jump-center right inaccessible region boundary #---> this point forward is universal for above inputs plot(c(a$theta1, a$theta1[1]), c(a$theta2, a$theta2[1]), type = "l", axes = FALSE, xlab = expression(theta), ylab = expression(kappa), main = "") # ensure axis labels are accurate points(a$theta1, a$theta2, lwd = 0.65) axis(side = 1, at = c(a$theta1hat, range(a$theta1))) axis(side = 2, at = c(a$theta2hat, range(a$theta2))) for (i in 1:nrow(pL)) { # ID & plot jump-center repair points in each quad m1R <- abs(a$theta1 - pR[i, 1]) / diff(range(a$theta1)) # match theta1 values m2R <- abs(a$theta2 - pR[i, 2]) / diff(range(a$theta2)) # match theta2 values m1L <- abs(a$theta1 - pL[i, 1]) / diff(range(a$theta1)) # match theta1 values m2L <- abs(a$theta2 - pL[i, 2]) / diff(range(a$theta2)) # match theta2 values fm <- which.min(m1R + m2R) + 1 # "right" points (nearest right boundary point) to <- which.min(m1L + m2L) - 1 # "to" points (nearest left boundary point) points(a$theta1[fm:to], a$theta2[fm:to], pch = 16, cex = 0.4, col = "blue") } points(pL, pch = 24, bg = "green") # points on left-side of inaccessible region points(pR, pch = 24, bg = "yellow") # points on right-side of inaccessible region points(pJ, pch = 24, bg = "red") # jump-center legend("topright", legend = c("CR boundary points", "jump-center (JC)", "JC left boundary", "JC right boundary", "JC repair points"), pch = c(1, rep(24, 3), 21), pt.bg = c("black", "red", "green", "yellow", "blue"), pt.cex = c(1, 1, 1, 1, 0.7), bty = "n", cex = 0.8) # next, record x (it will be used in the next chunk) and also complete the second plot x <- NULL x$theta <- a$theta1 x$kappa <- a$theta2 x$thetahat <- a$theta1hat # needed for next chunk x$kappahat <- a$theta2hat # needed for next chunk x$repair$q2jumpxy <- pJ[1, ] # needed for next chunk x$repair$q2jumpL <- pL[1, ] # needed for next chunk x$repair$q2jumpR <- pR[1, ] # needed for next chunk plot(c(x$theta, x$theta[1]), c(x$kappa, x$kappa[1]), type = "l", axes = FALSE, xlab = expression(theta), ylab = expression(kappa)) axis(side = 1, at = range(x$theta)) axis(side = 2, at = range(x$kappa)) ## ----fig.show = 'hold', warning = FALSE--------------------------------------- par(mar = c(4, 4, 1.5, 1)) # top-left confidence region extreme (zoomed-in) plot(x$theta, x$kappa, xlim = c(min(x$theta), x$thetahat), ylim = c(x$kappahat, max(x$kappa)), type = 'l', xlab = expression(theta), ylab = expression(kappa), main = "top-left (zoom-in)") # augment the top-left region graphic with plot points and jump-center information points(x$theta, x$kappa) jumppoints <- t(matrix(c(x$repair$q2jumpxy, x$repair$q2jumpL, x$repair$q2jumpR), nrow = 2)) points(jumppoints, pch = 24, bg = rep(c("red", "green", "yellow"))) # bottom-right confidence region extreme (zoomed-in) plot(x$theta, x$kappa, ylim = c(min(x$kappa), 0.95), type = 'l', xlab = expression(theta), ylab = expression(kappa), main = "bottom-right (zoom-in)") points(x$theta, x$kappa) ## ----fig.show = 'hold', warning = FALSE, eval = FALSE------------------------- # # top-left confidence region extreme (zoomed-in) # x <- crplot(bc, alpha = 0.1, distn = "gamma", cen = cen, showjump = TRUE, jumpinfo = TRUE, # jumpshift = c(0.5, 0.5, 0.5, 0.07), jumpuphill = c(0.01, 0.25, 0.01, 0.05), # xlim = c(min(x$theta), x$thetahat), ylim = c(x$kappahat, 4.15), # main = "top-left (zoom-in)") # #> [1] "Confidence region plot complete; made using 299 boundary points." # # # bottom-right confidence region extreme (zoomed-in) # plot(x$theta, x$kappa, ylim = c(min(x$kappa), 0.95), type = 'l', # xlab = expression(theta), ylab = expression(kappa), main = "bottom-right (zoom-in)") # points(x$theta, x$kappa) ## ----fig.show = 'hold', echo = FALSE------------------------------------------ par(mar = c(4, 4.5, 2, 1.5)) # results of above chunk are input directly here to save processing time on vignette (otherwise prohibitively long for CRAN) a$theta1 <- c(11092.63, 10433.13, 9807.861, 8634.626, 8073.472, 7796.363, 7520.319, 6952.024, 6371.019, 5763.591, 5061.823, 4510.083, 4126.239, 4077.722, 4069.232, 4065.99, 4065.685, 4065.64, 4065.629, 4065.625, 3256.963, 2686.811, 2189.392, 1995.414, 1795.015, 1425.904, 1337.262, 1306.203, 1299.947, 1297.033, 1294.258, 1293.585, 1292.921, 1292.265, 1291.618, 1290.979, 1290.662, 1290.348, 1290.035, 1289.725, 1289.416, 1289.263, 1289.11, 1288.958, 1288.806, 1288.654, 1288.503, 1288.353, 1288.203, 1288.054, 1287.905, 1287.756, 1287.608, 1287.461, 1287.314, 1287.167, 1287.021, 1286.876, 1286.731, 1286.586, 1286.442, 1286.156, 1285.871, 1285.588, 1285.308, 1285.029, 1284.752, 1284.614, 1284.477, 1283.939, 1283.41, 1282.887, 1282.372, 1281.362, 1280.381, 1276.731, 1270.661, 1262.765, 1259.821, 1373.132, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.422, 1557.423, 1557.423, 1557.424, 1557.425, 1557.426, 1557.429, 1557.433, 1557.437, 1557.444, 1557.45, 1557.469, 1557.489, 1557.527, 1557.603, 1557.723, 1557.844, 1558.102, 1558.36, 1558.877, 1559.69, 1560.503, 1561.782, 1563.064, 1565.083, 1567.108, 1578.992, 1591.141, 1616.27, 1653.549, 1693.321, 1753.816, 1820.462, 2097.608, 2518.007, 2826.445, 3253.11, 3937.509, 4439.756, 5098.404, 5620.529, 6212.872, 6854.903, 7520.319, 8305.757, 9067.957, 9792.367, 10483.48, 11730.83, 12880.77, 14784.04, 16510.45, 19141.57, 21541.62, 25392.41, 27192.46, 28930.35, 35822.83, 42440.97, 48878.38, 49595.8, 50309.54, 51411.59, 51837.43, 52262.06, 52591.09, 52919.42, 53173.97, 53428.11, 53625.22, 53822.09, 53974.84, 54127.44, 54245.87, 54364.21, 54456.07, 54547.88, 54619.15, 54690.39, 54745.7, 54801, 54843.93, 54865.39, 54886.85, 54894.62, 54902.38, 54908.34, 54914.29, 54918.86, 54923.43, 54926.93, 54930.43, 54933.12, 54935.81, 54937.87, 54939.93, 54941.51, 54943.09, 54944.31, 54945.52, 57217.34, 63771.85, 72528.15, 82785.21, 94037.64, 106050.1, 118733.2, 132062, 146038.9, 160677.3, 175994.4, 192006.8, 208727.9, 226164.7, 244313.4, 263153.7, 282636.9, 302665.4, 323049.4, 343404.4, 362852.3, 371710.1, 375523.5, 378740.9, 380484.8, 381161.3, 381674.1, 381725.1, 381772.9, 381858.6, 381930.6, 381961.1, 381987.9, 382010.8, 382020.8, 382029.7, 382037.6, 382044.5, 382050.3, 382055.1, 382057.1, 382058.8, 382060.2, 382061.3, 382062.2, 382062.7, 382063, 382063, 382062.7, 382062.1, 382061.2, 382060, 382058.6, 382056.8, 382052.3, 382046.5, 382039.5, 382031.2, 382021.6, 382010.7, 381998.4, 381969.6, 381935.1, 381894.6, 381795, 381668.6, 381513.2, 381326.3, 381220.1, 381163.6, 381104.8, 379695.7, 377165.9, 281718.9, 224935.4, 181174.9, 158666.2, 143303.7, 131713.7, 105104.6, 83459.78, 78636.93, 73979.37, 65086.54, 60833.62, 56668.29, 52572.15, 48513.22, 46299.82, 44078.33, 38943.62, 35777.59, 32403.98, 29268.35, 25363.77, 23505.32, 23144.91, 23109.1, 23100.73, 23098.79, 23098.34, 23098.17, 23098.14, 23098.12, 23098.12, 21951.15, 21922.4, 21223.25, 18865.36, 17491.4, 16398.69, 15484.77, 13878.77, 12551.42) a$theta2 <- c(2.142689, 2.180192, 2.218991, 2.302236, 2.347989, 2.372294, 2.397765, 2.454672, 2.520183, 2.598519, 2.705241, 2.805313, 2.885938, 2.896894, 2.898831, 2.899572, 2.899642, 2.899652, 2.899655, 2.899655, 3.115974, 3.320484, 3.554898, 3.666119, 3.795549, 4.0691, 4.128645, 4.142534, 4.144245, 4.144853, 4.145301, 4.145388, 4.145466, 4.145534, 4.145593, 4.145642, 4.145663, 4.145681, 4.145698, 4.145712, 4.145724, 4.145729, 4.145733, 4.145737, 4.14574, 4.145743, 4.145745, 4.145747, 4.145748, 4.145748, 4.145748, 4.145747, 4.145746, 4.145744, 4.145742, 4.145739, 4.145735, 4.145731, 4.145726, 4.145721, 4.145715, 4.145702, 4.145687, 4.145669, 4.14565, 4.145628, 4.145604, 4.145591, 4.145577, 4.145519, 4.145453, 4.145378, 4.145295, 4.145104, 4.144882, 4.143683, 4.139903, 4.12763, 4.089194, 3.786755, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502058, 3.502057, 3.502056, 3.502055, 3.502054, 3.502052, 3.502049, 3.502043, 3.502038, 3.502029, 3.502021, 3.501996, 3.501971, 3.501922, 3.501822, 3.501666, 3.501509, 3.501173, 3.500838, 3.500166, 3.499112, 3.498058, 3.496403, 3.494747, 3.492148, 3.489547, 3.474448, 3.459287, 3.428771, 3.385445, 3.341559, 3.278935, 3.215022, 2.993232, 2.747052, 2.609885, 2.459118, 2.278387, 2.176874, 2.069931, 2.000331, 1.933303, 1.871587, 1.816841, 1.761516, 1.715315, 1.676807, 1.644061, 1.592815, 1.552607, 1.497032, 1.455406, 1.40337, 1.364588, 1.314333, 1.294595, 1.277305, 1.221507, 1.181113, 1.149826, 1.146713, 1.14368, 1.139119, 1.137395, 1.135696, 1.134393, 1.133105, 1.132114, 1.131132, 1.130375, 1.129623, 1.129042, 1.128464, 1.128018, 1.127573, 1.127228, 1.126885, 1.126619, 1.126353, 1.126147, 1.125942, 1.125783, 1.125703, 1.125624, 1.125595, 1.125566, 1.125544, 1.125522, 1.125505, 1.125488, 1.125475, 1.125462, 1.125453, 1.125443, 1.125435, 1.125427, 1.125421, 1.125416, 1.125411, 1.125407, 1.117252, 1.096153, 1.072437, 1.049464, 1.028608, 1.010034, 0.993517, 0.9787697, 0.9655359, 0.9536073, 0.942821, 0.9330526, 0.9242103, 0.9162307, 0.909078, 0.9027452, 0.8972609, 0.892705, 0.8892422, 0.8872046, 0.8873349, 0.8887878, 0.8900543, 0.8918315, 0.8934649, 0.8944449, 0.8955606, 0.8957111, 0.8958643, 0.8961787, 0.8965045, 0.8966719, 0.8968423, 0.897016, 0.897104, 0.8971929, 0.8972826, 0.8973731, 0.8974645, 0.8975568, 0.8976033, 0.89765, 0.897697, 0.8977441, 0.8977915, 0.8978391, 0.897887, 0.8979351, 0.8979834, 0.898032, 0.8980808, 0.8981299, 0.8981792, 0.8982287, 0.8983286, 0.8984295, 0.8985314, 0.8986344, 0.8987384, 0.8988436, 0.8989498, 0.8991657, 0.8993862, 0.8996115, 0.9000775, 0.9005651, 0.9010762, 0.9016128, 0.9018913, 0.9020333, 0.9021772, 0.9049458, 0.9086121, 0.9878909, 1.040294, 1.09132, 1.123505, 1.148818, 1.170232, 1.229817, 1.29463, 1.31208, 1.330309, 1.369683, 1.391109, 1.414105, 1.439018, 1.466397, 1.482652, 1.50005, 1.545203, 1.577264, 1.61596, 1.657195, 1.717887, 1.751489, 1.758431, 1.759129, 1.759292, 1.75933, 1.759339, 1.759342, 1.759343, 1.759343, 1.759343, 1.782534, 1.783137, 1.79814, 1.854264, 1.891683, 1.924498, 1.954384, 2.013396, 2.069931) a$theta1hat <- 7520.318946 # MLE a$theta2hat <- 2.069931 # MLE pJ <- matrix(c(2189.392, 3.115974, 48513.22, 1.170232), nrow = 2, byrow = TRUE) # jump-center location pL <- matrix(c(1557.422, 3.502058, 21922.4, 1.783137), nrow = 2, byrow = TRUE) # jump-center left inaccessible region boundary pR <- matrix(c(4065.629, 2.899655, 54943.09, 1.125416), nrow = 2, byrow = TRUE) # jump-center right inaccessible region boundary #---> this point forward is universal for above inputs plot(c(a$theta1, a$theta1[1]), c(a$theta2, a$theta2[1]), type = "l", axes = FALSE, xlim = c(min(x$theta), x$thetahat), ylim = c(x$kappahat, 4.15), xlab = expression(theta), ylab = expression(kappa), main = "top-left (zoom-in)") # ensure axis labels are accurate points(a$theta1, a$theta2, lwd = 0.65) axis(side = 1, at = c(a$theta1hat, min(a$theta1))) axis(side = 2, at = c(a$theta2hat, max(a$theta2))) for (i in 1:nrow(pL)) { # ID & plot jump-center repair points in each quad m1R <- abs(a$theta1 - pR[i, 1]) / diff(range(a$theta1)) # match theta1 values m2R <- abs(a$theta2 - pR[i, 2]) / diff(range(a$theta2)) # match theta2 values m1L <- abs(a$theta1 - pL[i, 1]) / diff(range(a$theta1)) # match theta1 values m2L <- abs(a$theta2 - pL[i, 2]) / diff(range(a$theta2)) # match theta2 values fm <- which.min(m1R + m2R) + 1 # "right" points (nearest right boundary point) to <- which.min(m1L + m2L) - 1 # "to" points (nearest left boundary point) points(a$theta1[fm:to], a$theta2[fm:to], pch = 16, cex = 0.4, col = "blue") } points(pL, pch = 24, bg = "green") # points on left-side of inaccessible region points(pR, pch = 24, bg = "yellow") # points on right-side of inaccessible region points(pJ, pch = 24, bg = "red") # jump-center legend("topright", legend = c("CR boundary points", "jump-center (JC)", "JC left boundary", "JC right boundary", "JC repair points"), pch = c(1, rep(24, 3), 21), pt.bg = c("black", "red", "green", "yellow", "blue"), pt.cex = c(1, 1, 1, 1, 0.7), bty = "n", cex = 0.8) points(a$theta1hat, a$theta2hat, pch = 3) # complete next plot: x <- NULL x$theta <- a$theta1 x$kappa <- a$theta2 x$thetahat <- a$theta1hat # needed for next chunk x$kappahat <- a$theta2hat # needed for next chunk # bottom-right confidence region extreme (zoomed-in) plot(x$theta, x$kappa, ylim = c(min(x$kappa), 0.95), type = 'l', xlab = expression(theta), ylab = expression(kappa), main = "bottom-right (zoom-in)") points(x$theta, x$kappa) ## ----fig.width = 6.2, fig.height = 6, fig.show = 'hold'----------------------- # final plot, with all necessary repairs complete plot(c(x$theta, x$theta[1]), c(x$kappa, x$kappa[1]), type = "l", axes = FALSE, xlab = expression(theta), ylab = expression(kappa), main = "gamma 90% CR for BF data (final)") axis(side = 1, at = range(x$theta)) axis(side = 2, at = range(x$kappa)) points(x$thetahat, x$kappahat, pch = 3) ## ----------------------------------------------------------------------------- #> [1] "--------------------------------------------------------------------------------------" #> [1] "R uniroot failure searching for confidence region boundary---challenging parameters or shape." #> [1] "Unable to produce a confidence region for the given sample and/or parameterization. " #> [1] "--------------------------------------------------------------------------------------" ## ----fig.show = 'hold', warning = FALSE, error = FALSE, eval = FALSE---------- # # crplot is unable to plot this 98% confidence region # crplot(dataset = c(1.5, 2), alpha = 0.01, distn = "invgauss") # #> [1] "--------------------------------------------------------------------------------------" # #> [1] "R uniroot failure searching for confidence region boundary---challenging parameters and/or shape." # #> [1] "Unable to produce a confidence region for the given sample and/or parameterization. " # #> [1] "--------------------------------------------------------------------------------------" ## ----fig.show = 'hold', warning = FALSE--------------------------------------- # a plot without jump-center repairs is attainable, but its 3rd and 4th quadrants, # relative to the MLE, are in need of jump-center repairs crplot(dataset = c(1.5, 2), alpha = 0.01, distn = "invgauss", repair = FALSE, sf = c(3, 3), main = "without repairs") # a complete plot is returned by increasing alpha to values >= 0.03 crplot(dataset = c(1.5, 2), alpha = 0.05, distn = "invgauss", main = "95% CR", sf = c(3, 3)) # adjusting jumpshift and jumpuphill parameters x <- crplot(dataset = c(1.5, 2), alpha = 0.02, "invgauss", jumpinfo = TRUE, sf = c(3, 3), jumpshift = 0.1, jumpuphill = 0.001, main = "98% CR") plot(c(x$mu, x$mu[1]), c(x$lambda, x$lambda[1]), type = "l", axes = FALSE, xlab = expression(mu), ylab = expression(lambda), main = "98% CR") axis(side = 1, at = format(range(x$mu), digits = 3)) axis(side = 2, at = format(range(x$lambda), digits = 3))