You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
#' Estimate Pareto Parameters#'#' @family Parameter Estimation#' @family Pareto#'#' @author Steven P. Sanderson II, MPH#'#' @details This function will attempt to estimate the Pareto shape and scale#' parameters given some vector of values.#'#' @description The function will return a list output by default, and if the parameter#' `.auto_gen_empirical` is set to `TRUE` then the empirical data given to the#' parameter `.x` will be run through the `tidy_empirical()` function and combined#' with the estimated Pareto data.#'#' Two different methods of shape parameters are supplied:#' - LSE#' - MLE#'#' @param .x The vector of data to be passed to the function.#' @param .auto_gen_empirical This is a boolean value of TRUE/FALSE with default#' set to TRUE. This will automatically create the `tidy_empirical()` output#' for the `.x` parameter and use the `tidy_combine_distributions()`. The user#' can then plot out the data using `$combined_data_tbl` from the function output.#'#' @examples#' library(dplyr)#' library(ggplot2)#'#' x <- mtcars[["mpg"]]#' output <- util_pareto1_param_estimate(x)#'#' output$parameter_tbl#'#' output$combined_data_tbl |>#' tidy_combined_autoplot()#'#' set.seed(123)#' t <- tidy_pareto1(.n = 100, .shape = 1.5, .min = 1)[["y"]]#' util_pareto1_param_estimate(t)$parameter_tbl#'#' @return#' A tibble/list#'#' @name util_pareto1_param_estimateNULL#' @export#' @rdname util_pareto1_param_estimateutil_pareto1_param_estimate<-function(.x, .auto_gen_empirical=TRUE) {
# Tidyeval ----x_term<- as.numeric(.x)
minx<- min(x_term)
maxx<- max(x_term)
n<- length(x_term)
unique_terms<- length(unique(x_term))
# Checks ----if (!is.vector(x_term, mode="numeric") || is.factor(x_term)) {
rlang::abort(
message="'.x' must be a numeric vector.",
use_cli_format=TRUE
)
}
if (n<2|| any(x_term<=0) ||unique_terms<2) {
rlang::abort(
message="'.x' must contain at least two non-missing distinct values. All values of '.x' must be positive.",
use_cli_format=TRUE
)
}
# Get params ----# LSEppc<-0.375fhat<-stats::ppoints(n, a=ppc)
lse_coef<-stats::lm(log(1-fhat) ~ log(sort(x_term)))$coefficientslse_shape<--lse_coef[[2]]
lse_min<- exp(lse_coef[[1]] /lse_shape)
# MLEmle_min<- min(x_term)
mle_shape<-n/ sum(log(x_term/mle_min))
# Return Tibble ----if (.auto_gen_empirical) {
te<- tidy_empirical(.x=x_term)
td_lse<- tidy_pareto1(.n=n, .shape= round(lse_shape, 3), .min= round(lse_min, 3))
td_mle<- tidy_pareto1(.n=n, .shape= round(mle_shape, 3), .min= round(mle_min, 3))
combined_tbl<- tidy_combine_distributions(te, td_lse, td_mle)
}
ret<-dplyr::tibble(
dist_type= rep("Pareto", 2),
samp_size= rep(n, 2),
min= rep(minx, 2),
max= rep(maxx, 2),
method= c("LSE", "MLE"),
est_shape= c(lse_shape, mle_shape),
est_min= c(lse_min, mle_min)
)
# Return ----
attr(ret, "tibble_type") <-"parameter_estimation"
attr(ret, "family") <-"pareto"
attr(ret, "x_term") <-.x
attr(ret, "n") <-nif (.auto_gen_empirical) {
output<-list(
combined_data_tbl=combined_tbl,
parameter_tbl=ret
)
} else {
output<-list(
parameter_tbl=ret
)
}
return(output)
}
#' Calculate Akaike Information Criterion (AIC) for Pareto Distribution#'#' This function calculates the Akaike Information Criterion (AIC) for a Pareto distribution fitted to the provided data.#'#' @family Utility#' @family Pareto#' @author Steven P. Sanderson II, MPH#'#' @description#' This function estimates the shape and scale parameters of a Pareto distribution#' from the provided data using maximum likelihood estimation,#' and then calculates the AIC value based on the fitted distribution.#'#' @param .x A numeric vector containing the data to be fitted to a Pareto distribution.#'#' @details#' This function fits a Pareto distribution to the provided data using maximum#' likelihood estimation. It estimates the shape and scale parameters#' of the Pareto distribution using maximum likelihood estimation. Then, it#' calculates the AIC value based on the fitted distribution.#'#' Initial parameter estimates: The function uses the method of moments estimates#' as starting points for the shape and scale parameters of the Pareto distribution.#'#' Optimization method: The function uses the optim function for optimization.#' You might explore different optimization methods within optim for potentially#' better performance.#'#' Goodness-of-fit: While AIC is a useful metric for model comparison, it's#' recommended to also assess the goodness-of-fit of the chosen model using#' visualization and other statistical tests.#'#' @examples#' # Example 1: Calculate AIC for a sample dataset#' set.seed(123)#' x <- tidy_pareto1()$y#' util_pareto_aic(x)#'#' @return#' The AIC value calculated based on the fitted Pareto distribution to the provided data.#'#' @name util_pareto_aicNULL#' @export#' @rdname util_pareto_aicutil_pareto_aic<-function(.x) {
# Tidyevalx<- as.numeric(.x)
n<- length(x)
# Negative log-likelihood function for Pareto distributionneg_log_lik_pareto<-function(par, data) {
shape<-par[1]
min<-par[2]
-sum(actuar::dpareto1(data, shape=shape, min=min, log=TRUE))
}
# Get initial parameter estimates: method of momentspe<-TidyDensity::util_pareto1_param_estimate(x)$parameter_tbl|>
subset(method=="MLE")
# Fit Pareto distribution using optimfit_pareto<-stats::optim(
c(pe$est_shape, pe$est_min),
neg_log_lik_pareto,
data=x
)
# Extract log-likelihood and number of parameterslogLik_pareto<--fit_pareto$valuek_pareto<-2# Number of parameters for Pareto distribution (shape and min)# Calculate AICAIC_pareto<-2*k_pareto-2*logLik_pareto# Return AICreturn(AIC_pareto)
}
#' Distribution Statistics for Pareto1 Distribution#'#' @family Pareto#' @family Distribution Statistics#'#' @details This function will take in a tibble and returns the statistics#' of the given type of `tidy_` distribution. It is required that data be#' passed from a `tidy_` distribution function.#'#' @description Returns distribution statistics in a tibble.#'#' @param .data The data being passed from a `tidy_` distribution function.#'#' @examples#' library(dplyr)#'#' tidy_pareto1() |>#' util_pareto1_stats_tbl() |>#' glimpse()#'#' @return#' A tibble#'#' @name util_pareto1_stats_tblNULL#' @export#' @rdname util_pareto1_stats_tblutil_pareto1_stats_tbl<-function(.data) {
# Immediate check for tidy_ distribution functionif (!"tibble_type"%in% names(attributes(.data))) {
rlang::abort(
message="You must pass data from the 'tidy_dist' function.",
use_cli_format=TRUE
)
}
if (attributes(.data)$tibble_type!="tidy_pareto_single_parameter") {
rlang::abort(
message="You must use 'tidy_pareto1()'",
use_cli_format=TRUE
)
}
# Datadata_tbl<-dplyr::as_tibble(.data)
atb<- attributes(data_tbl)
xm<-atb$.minalpha<-atb$.shapestat_mean<- ifelse(alpha<=1, Inf, (alpha*xm) / (alpha-1))
stat_mode<-xmstat_coef_var<- ifelse(
alpha<=2,
Inf,
sqrt((alpha) / ((alpha-1)^2* (alpha-2)))
)
stat_sd<- ifelse(
alpha<=1,
Inf,
sqrt((alpha*xm^2) / ((alpha-1)^2* (alpha-2)))
)
stat_skewness<- ifelse(
alpha<=3,
"undefined",
(2* (1+alpha)) / (alpha-3) * sqrt((alpha-2) /alpha)
)
stat_kurtosis<- ifelse(
alpha<=4,
"undefined",
(6* (alpha^3+alpha^2-6*alpha-2)) / (alpha* (alpha-3) * (alpha-4))
)
# Data Tibbleret<-dplyr::tibble(
tidy_function=atb$tibble_type,
function_call=atb$dist_with_params,
distribution="Pareto1",
distribution_type="Continuous",
points=atb$.n,
simulations=atb$.num_sims,
mean=stat_mean,
mode_lower=stat_mode,
range= paste0(xm, " to Inf"),
std_dv=stat_sd,
coeff_var=stat_coef_var,
skewness=stat_skewness,
kurtosis=stat_kurtosis,
computed_std_skew= tidy_skewness_vec(data_tbl$y),
computed_std_kurt= tidy_kurtosis_vec(data_tbl$y),
ci_lo= ci_lo(data_tbl$y),
ci_hi= ci_hi(data_tbl$y)
)
# Returnreturn(ret)
}
Param Estimates
Function:
Example:
AIC Function
Function:
Example:
Stats Tibble
Function:
Example:
The text was updated successfully, but these errors were encountered: