--- title: "crplot Advanced Options" author: "Christopher Weld , Lawrence Leemis " date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{crplot_advanced} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ## Introduction This is the second of two vignettes on the `crplot` function. It details advanced features which are helpful when troubleshooting plot issues, focusing on its jump-center algorithm parameters. For a basic comprehensive summary of `crplot` usage, please reference the first vignette titled *crplot* (users should familiarize with it prior to viewing this vignette). Both vignettes are available via links on the [conf](https://CRAN.R-project.org/package=conf) package webpage. The `crplot` function generates a two-dimensional confidence region plot for the specified two-parameter parametric distribution, fitted to a dataset. Details of the `crplot` algorithm are available in its corresponding publication^[Weld, C., Loh, A., and Leemis, L. (in press), "Plotting Likelihood-Ratio Based Confidence Regions for Two-Parameter Univariate Probability Models", *The American Statistician*]. ## Jump-Center Components A jump-center is necessary to plot confidence region boundary points that are inaccessible via the radial profile technique (plotting algorithm of `crplot`) because a line joining them to the MLE crosses the confidence region boundary multiple times. The jump-center algorithm "repairs" the inaccessible region by providing an alternate point-of-reference (or alternate "center-point", hence its name) from which to execute the radial profile plotting technique. See `crplot` vignette 1 of 2 for additional details. Many advanced `crplot` options involve customizing its jump-center parameters. Before exploring those methods, it is first important to familiarize with jump-center components. The plot annotations below will aid in those definitions. ```{r, 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) ``` *(note: select word coloring below is chosen to correspond with the appropriate reference points above)* * *Jump-Center \& its Left and Right Boundaries*. Customization is informed by the default jump-center (JC) and its JC left and JC right boundary reference points (pictured above). The "left" and "right" directional labels represent their location relative to the jump-center (when looking "out" from the jump-center into the inaccessible region). Maximum angle tolerance constraints are not satisfied at the left and/or right inaccessible region border points, and the jump-center repair aims to plot points within this uncharted region. * *Quadrant*. The jump-center position with respect to the MLE ($+$) is significant. As a convenience, quadrant notation relative to the MLE is used. For example, the jump-center pictured above lies in the fourth quadrant relative to the MLE. *All future references to quadrants within this vignette adopt this perspective (quadrants defined as if the MLE were the origin)*. * *Gap*. Two inaccessible region orientations are possible: a vertical (shown above) or horizontal gap. The vertical gap shown above has an inaccessible region orientated counterclockwise from the MLE ($+$). For a jump-center in the fourth quadrant, a clockwise inaccessible region orientation corresponds to a horizontally oriented gap. * *Jump-Shift*. The direction of the jump-center, relative to the MLE ($+$), intercepts the horizontal or vertical gap according to a jump-shift parameter. The jump-shift parameter represents the fraction of its gap where an azimuth to the jump-center crosses. The default, pictured above, is a jump-shift = 0.5 (or bisecting its gap). For this example, decreasing this value moves the azimuth from the MLE to the jump-center nearer its right boundary reference point. Increasing its value would move it nearer the left boundary. * *Jump-Uphill*. The jump-center must locate *inside* of the existing confidence region boundary. Considering the log-likelihood function governs the landscape from which the confidence region contour is found, points inside its region are those "uphill" from its border. The jump-uphill parameter quantifies how far uphill (given as a measure of level of significance) from its confidence region boundary the jump-center locates, and thus determines its distance from the MLE (or confidence region boundary). For example, consider an $89\%$ versus a $90\%$ confidence region (or levels of significance of $0.11$ and $0.1$, respectively). The $89\%$ confidence region is smaller; its boundary is inside that of the $90\%$ region and represents a "jump-uphill" measure of $0.01$ (in level of significance). The default parameterization for jump-uphill is $\min(\alpha, 0.01)$. ## Plot errors Plot issues using `crplot` are relatively uncommon, but possible. Certain datasets and/or distributions pose greater plot challenges than others. In particular, small datasets and/or heavily censored datasets can result in unusual (strong deviations from elliptical, even to the point of non-convex) shapes, making their boundary difficult to map using the radial profile log likelihood technique^[Jaeger, A. (2016), "Computation of Two- and Three-Dimensional Confidence Regions With the Likelihood Ratio", *The American Statistician*, 70, 395--398] employed by `crplot`. When the `crplot` algorithm is stretched beyond its default limitations, adjusting a series of optional arguments can further extend its practical reach. Optional arguments in the `crplot` function enable algorithm customization, which is often enough to overcome plot issues when they arise. ## Example This heavily right-censored example produces a particularly challenging confidence region shape, enabling us to illustrate a suite of optional arguments available to its user to customize the `crplot` algorithm. Understanding the impact of optional argument adjustments will facilitate their use elsewhere when appropriate. #### Heavily Right-Censored Dataset Bearing cage fracture data is taken from Abernethy et. al. (1983)^[Abernethy, R. B., Breneman, J. E., Medlin, C. H., and Reinman, G. L. (1983), "Weibull Analysis Handbook", *Aero Propulsion Laboratory*, 43--47]. Among 1703 samples there are six observed failures. The remaining 1697 right-censored samples occur over a variety of right-censored values between 50--2050 hours. Both data subsets are given below with the variables `bc_obs` and `bc_cen`, respectively. ```{r} 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)) ``` #### Default Confidence Region Results Confidence regions corresponding to a variety of parametric distributions for this heavily right-censored dataset are given below. Each plot uses a level of significance of `alpha = 0.1` (a 90\% confidence region). ```{r, 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])) } ``` ```{r, 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) } } ``` #### Investigating Plot Issues The gamma distribution confidence region plot for the bearing cage fracture dataset takes a suspect shape. Nearer the bottom of its plot there are long segments that are seemingly missing confidence region boundary points. Relatively sharp vertex angles flank each of these line segments, a characteristic indicative of a plot issue. To assess what is happening, jump-center reference points are plot using `showjump = TRUE`. Additionally, its plot information is returned using the optional argument `jumpinfo = TRUE`. The `jumpinfo` argument is used in lieu of its alternative `info = TRUE` because we want jump-center repair information returned in addition to the final confidence region plot information. It is then available for subsequent use, as shown by plotting those points on the "without repairs" plot below. ```{r, 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" ``` ```{r, 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") ``` Setting `jumpinfo = TRUE` returns a `repair` attribute, which is a list relevant to all jump-center repairs. The `repair` labels are output to the console above, and given in greater detail next. Prefixes associated with `repair` are: - `q1`, `q2`, `q3`, `q4`. Each label begins with an abbreviation specifying what quadrant (with respect to the MLE) the jump-center repair is made. For example, a `q1` prefix indicates a jump-center repair located in the 1st quadrant (the jump-center horizontal and vertical values are both greater than the MLE's respective values). In our example, `q2` and `q4` prefixes indicate jump-center repairs in the 2nd and 4th quadrant relative to the MLE (to its upper-left and lower-right, respectively). Suffixes associated with `repair` are: - `jumpuphill` quantifies the amount uphill (with regards to level of significance) from the confidence region boundary that the jump-center will locate, - `jumpshift` quantifies the point along the vertical or horizontal "gap" (as a fractional value between 0 and 1) that the azimuth from the MLE to the jump-center will cross, - `jumpxy` stores the coordinates for the jump-center, - `jumpL` stores the left (with perspective taken from the jump-center location towards its inaccessible region) confidence region boundary point bordering its inaccessible region, and - `jumpR` stores the left (with perspective taken from the jump-center location towards its inaccessible region) confidence region boundary point bordering its inaccessible region. Seeing its jump-center location enables a targeted strategy---using the optional arguments given next---to improve the plot when necessary. ### Adjusting Optional Arguments Optional arguments whose adjustment can refine `crplot` results include: - `jumpuphill` - `jumpshift` - `ellipse_n` - `maxcount` Optional argument adjustments can be applied separately or in combination, each serving a unique purpose. Combinations of `jumpshift` and/or `jumpuphill` are the recommended first option to address plot regions that remain inaccessible to its jump-center. The fourth quadrant repairs (with respect to the MLE) of the aforementioned example will illustrate optional argument uses below since they are more easily visible than its second quadrant repairs. #### $\circ$ `jumpuphill` Increasing this optional argument from its default value, `jumpuphill = min(alpha, 0.01)`, moves the jump-center nearer the MLE and further from its confidence region boundary. Doing so in our example will enable the fourth quadrant jump-center to achieve a better line-of-sight to the bottom portion of its confidence region (where a perceived lack of plot points appears problematic). ```{r, 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." ``` ```{r, 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) ``` The confidence region's bottom-most curve benefits substantially from this `jumpuphill` adjustment, but a problematic area (lacking plot points) remains in the upper curve at higher $\theta$ values. #### $\circ$ `jumpshift` Adjusting this optional argument from its default value, `jumpshift = 0.5`, shifts the radial azimuth (with respect to the MLE) to its jump-center. Decreasing its value generally "pushes" the jump-center location further along its confidence region boundary, nearer the inaccessible region. Doing so in our example enables the fourth quadrant jump-center to achieve a better line-of-sight to the upper portion of its confidence region curve at higher $\theta$ values (where a perceived lack of plot points appears problematic). ```{r, 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." ``` ```{r, 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) ``` The confidence region's upper curve (at higher $\theta$ values) benefits substantially from this `jumpshift` adjustment, but does so at the expense of its bottom curve. The bottom curve (at minimal $\kappa$ values) plot difficulties are compounded by the perspective taken by this new jump-center location. *Caution*: although decreasing `jumpshift` generally improves access to inaccessible regions, inevitably a `jumpshift` value too near zero will result in multiple solutions for its jump-center at the specified `jumpuphill` value. Under those circumstance, the point nearer its MLE is typically returned, resulting in a jump-center whose location is not between the inaccessible boundary points, and is therefore unable to plot its inaccessible region. Using `jumpshift = 0.01` rather than `jumpshift = 0.1` with the syntax above demonstrates this undesirable paradigm, and is left as an exercise to the reader. Another example of a `jumpshift` value that is prohibitively low for its circumstance is given later in this vignette, in the *Optional Argument Combinations* section. #### $\circ$ `ellipse_n` Adjusting this optional argument from its default value, `ellipse_n = 4`, combines the elliptically-oriented plot technique with the default smoothing-search algorithm. Details of the elliptically-oriented plot algorithm are given in the first `crplot` vignette. It increases the final number of plot points, albeit at a computational expense. Used here, it improves but does not entirely correct our example's plot concerns. ```{r, 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." ``` ```{r, 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) ``` #### $\circ$ `maxcount` The `maxcount` parameter specifies the number of smoothing-search iterations `crplot` performs before terminating its algorithm. The default value (30) is typically sufficient to plot all accessible regions to its maximum degree tolerance constraint (`maxdeg`), so reaching this limit is indicative of an inaccessible region plot issue. On rare occasion however, increasing `maxcount` may improvement select plot regions, albeit at a noteworthy computational cost. This is true for our example plot, whose bottom region continues to (slowly) populate with relevant points as `maxcount` increases. ```{r, 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." ``` ```{r, 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))) ``` Increasing `maxcount` is effective at identifying additional relevant points here due to poor circumstance surrounding its fourth quadrant jump-center position. The smoothing search algorithm iteratively locates additional plot points based on azimuths (or angles) associated with the mid-point of its existing points (wherever `maxdeg` angle constraints are not yet met). The mid-point between adjacent points along its bottom curve are so offset and spread from the jump-center that, despite aiming towards the mid-point of its long vacant region, additional points are successively added much nearer its very far-right existing point. #### $\circ$ Optional argument combinations Combinations of the aforementioned optional arguments are typically capable of correcting the most challenging plot circumstance. For the bearing cage fracture data example, `jumpuphill` and `jumpshift` parameter adjustments---influencing the jump-center distance and direction from the MLE, respectively---prove sufficient to adequately repair plot issues in its fourth quadrant. One solution, using `jumpshift = 0.07` and `jumpuphill = 0.05`, is given below. ```{r, 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)) ``` ```{r, 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)) ``` A closer review of its lower-right and upper-left plot regions is appropriate to ensure those extremes are adequately captured. Those plots reveal that the quadrant four repairs were successful, but applying its customizations uniformly for all jump-center locations results in an inaccurate plot region for quadrant two. It encounters an issue discussed as "*Caution*: ..." in the `jumpshift` section above. Specifically, its decreased `jumpshift` value causes the MLE-to-jump-center azimuth to pass too near its inaccessible region right boundary. It subsequently fails to locate its jump-center by the "far" confidence region boundary, and instead crosses another solution for the `jumpuphill` parameter before its inaccessible region. ```{r, 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) ``` Different `jumpshift` and `jumpuphill` values for quadrant two and four repairs are necessary to adequately plot each of its extreme regions to a high level of detail. Itemizing customizations in this manner is possible, and detailed in the next section. #### Varying `jumpshift` and `jumpuphill` values by quadrant. A plot with multiple jump-center repairs may warrant a different `jumpshift` and `jumpuphill` value for each respective jump-center. Given scalar `jumpshift` and/or `jumpuphill` arguments, the `crplot` algorithm assumes its value applies regardless of the jump-center repair location (or quadrant). To customize these parameters to jump-center location-specific values, assign `jumpshift` and/or `jumpuphill` a vector of length four; its values then correspond uniquely to the four respective quadrants surrounding the MLE. When providing a quadrant-specific vector of `jumpshift` or `jumpuphill` values, any placeholder value is acceptable for quadrants with no jump-center. For our example, quadrants one and three contain no jump-center repairs and we can therefore arbitrarily assign the default values of 0.5 and 0.01 for `jumpshift` and `jumpuphill` in those respective locations. Through separate analysis (analogous to that for quadrant four), we determine that the quadrant two jump-center adequately plots its associated region using `jumpuphill = 0.25` (keeping the default value for `jumpshift = 0.5`). Notice how the maximum $\kappa$ value increases as a result. This solution is combined with the quadrant four solution (`jumpshift = 0.07` and `jumpuphill = 0.05`) below. ```{r, 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) ``` ```{r, 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) ``` Satisfied with plot accuracy in all regions, the final $90\%%$ confidence region for the heavily right-censored bearing fracture data is given: ```{r, 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) ``` ## Trouble-shooting R Uniroot Errors Albeit uncommon, if `crplot` is unable to return a confidence region plot the following message prints to the console: ```{r} #> [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] "--------------------------------------------------------------------------------------" ``` The R `uniroot` function numerically solves for confidence region boundary points within `crplot`. No plot returns when this numeric method fails. Failures can occur due-to, or despite `uniroot` parameters `upper` and `lower` properly bracketing the confidence region (uniroot solves using a bisection method). This typically occurs under the most challenging plot circumstance such as heavily censored and/or small sample sizes, or scales that approach the limits of R's finite precision arithmetic. Good options to troubleshoot this error include: increasing `alpha`, setting `repair = FALSE`, and adjusting the `jumpshift` and/or `jumpuphill` parameters. Increasing `alpha` is effective because it reduces the confidence region size, and this is can remedy `uniroot` upper-bound issues. The latter two are effective options because both influence the jump-center repairs algorithm. The jump-center repairs algorithm is usually the source of uniroot errors since it addresses the most "extreme" corners of the confidence region plot, where `uniroot` numeric issues tend to occur. Setting `repair = FALSE` may enable return of a working-copy of the confidence region (complete with exception of its radially inaccessible regions). To attain a complete confidence region, adjusting `jumpshift` and/or `jumpuphill` inherently will assign its jump-center a different location, from where its `uniroot` calculations may not experience similar numeric difficulties. The example below demonstrates circumstance resulting in a plot error, and attainable solutions using the strategy described above. ```{r, 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] "--------------------------------------------------------------------------------------" ``` The three troubleshooting options given above are illustrated next. For this example, each successfully produces a confidence region plot. ```{r, 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)) ```