XLSX with 300WSs added, INFO boxes added to OPTIMIZE
This commit is contained in:
1 parent
e263ba66e4
commit
be376b096c
3 files changed
+221
-100
No files matched your search
+44
@@ -1626,3 +1626,47 @@ bendFUN <- function(coeffsMR, coeffsMU) {
|
||||
)
|
||||
return(BPsMR_MU)
|
||||
}
|
||||
|
||||
#' Fitting of Gumbel location and scale parameters, and plotting function
|
||||
#'
|
||||
#'
|
||||
#' @param H Histogram density results with the mid-points of the histograms.
|
||||
#' @returns A ggplot object and a table of Gumbel quantiles 0.95, 0.99, and 0.999.
|
||||
#' @export
|
||||
#' @examples
|
||||
#' suppressMessages(source("../../dev/setup.R"))
|
||||
#'
|
||||
#' H = structure(list(breaks = c(1.02, 1.04, 1.06, 1.08, 1.1, 1.12, 1.14),
|
||||
#' counts = c(2L, 12L, 11L, 4L, 0L, 1L),
|
||||
#' density = c(3.33333333333333, 20, 18.3333333333333, 6.66666666666666, 0, 1.66666666666667),
|
||||
#' mids = c(1.03, 1.05, 1.07, 1.09, 1.11, 1.13), xname = 'UppAsympt',
|
||||
#' equidist = TRUE), class = "histogram")
|
||||
#' Titel <- "test"
|
||||
#'
|
||||
#' GumbPlotFUN(H)
|
||||
GumbPlotFUN <- function(H, Titel="Maximum deviation of confidence limit") {
|
||||
|
||||
H_df <- data.frame(Density = H[[3]], Mids=H[[4]])
|
||||
startlist <- list(mu = 1.1, b = 1)
|
||||
|
||||
GuM <- gsl_nls(
|
||||
fn = Density ~ 1/b*exp(-((Mids-mu)/b+exp(-(Mids-mu)/b))),
|
||||
data = H_df,
|
||||
start = startlist, # race=T,
|
||||
control = gsl_nls_control(xtol = 1e-6, ftol = 1e-6, gtol = 1e-6)
|
||||
)
|
||||
SuGuM <- summary(GuM)$coefficients
|
||||
X_gu <- seq(min(H[["breaks"]]), max(H[["breaks"]]), by=0.001)
|
||||
DGu <- dgumbel(X_gu, loc=SuGuM[1,1], scale = SuGuM[2,1])
|
||||
DGuDF <- data.frame(X = X_gu, Y = DGu)
|
||||
Quantiles <- qgumbel(c(0.95,0.99,0.999),loc=SuGuM[1,1], scale = SuGuM[2,1] )
|
||||
INV_Quantiles <- 1/Quantiles
|
||||
|
||||
p_gu <- ggplot(H_df) +
|
||||
geom_col(mapping=aes(x=Mids, y=Density)) +
|
||||
geom_line(DGuDF, mapping=aes(x=X, y=Y), inherit.aes = F) +
|
||||
geom_vline(xintercept = Quantiles, linetype=2, col="blue") +
|
||||
ggtitle(Titel) +
|
||||
theme_bw()
|
||||
return(list(p_gu, data.frame(Probabilities= c(0.95,0.99,0.999), INV_Quantiles,Quantiles)))
|
||||
}
|
||||
Reference in new issue
Block a user