Skip to content

Commit 4628d68

Browse files
authored
Merge pull request #438 from spsanderson/development
Fixes #437
2 parents a7e762e + 2538009 commit 4628d68

9 files changed

+421
-336
lines changed

NEWS.md

+3
Original file line numberDiff line numberDiff line change
@@ -9,12 +9,15 @@ None
99
3. Fix #414 - Add function `util_chisquare_param_estimate()` to estimate the parameters of the chi-square distribution.
1010
4. Fix #417 - Add function `tidy_mcmc_sampling()` to sample from a distribution using MCMC.
1111
This outputs the function sampled data and a diagnostic plot.
12+
5. Fix #420 - Add functions `util_dist_aic()` functions to calculate the AIC for a distribution.
1213

1314
## Minor Fixes and Improvements
1415
1. Fix #401 - Update `tidy_multi_single_dist()` to respect the `.return_tibble` parameter
1516
2. Fix #406 - Update `tidy_multi_single_dist()` to exclude the `.return_tibble` parameter
1617
from returning in the distribution parameters.
1718
3. Fix #413 - Update documentation to include `mcmc` where applicable.
19+
4. Fix #240 - Update `tidy_distribution_comparison()` to include the new AIC calculations
20+
from the dedicated `util_dist_aic()` functions.
1821

1922
# TidyDensity 1.3.0
2023

R/est-param-chisq.R

+43-28
Original file line numberDiff line numberDiff line change
@@ -95,36 +95,51 @@ util_chisquare_param_estimate <- function(.x, .auto_gen_empirical = TRUE) {
9595
}
9696

9797
# Parameters ----
98-
estimate_chisq_params <- function(data) {
99-
# Negative log-likelihood function
100-
negLogLik <- function(df, ncp) {
101-
-sum(stats::dchisq(data, df = df, ncp = ncp, log = TRUE))
102-
}
103-
104-
# Initial values (adjust based on your data if necessary)
105-
start_vals <- list(df = trunc(var(data)/2), ncp = trunc(mean(data)))
106-
107-
# MLE using bbmle
108-
mle_fit <- bbmle::mle2(negLogLik, start = start_vals)
109-
# Return estimated parameters as a named vector
110-
df <- dplyr::tibble(
111-
est_df = bbmle::coef(mle_fit)[1],
112-
est_ncp = bbmle::coef(mle_fit)[2]
113-
)
114-
return(df)
98+
# estimate_chisq_params <- function(data) {
99+
# # Negative log-likelihood function
100+
# negLogLik <- function(df, ncp) {
101+
# -sum(stats::dchisq(data, df = df, ncp = ncp, log = TRUE))
102+
# }
103+
#
104+
# # Initial values (adjust based on your data if necessary)
105+
# start_vals <- list(df = trunc(var(data)/2), ncp = trunc(mean(data)))
106+
#
107+
# # MLE using bbmle
108+
# mle_fit <- bbmle::mle2(negLogLik, start = start_vals)
109+
# # Return estimated parameters as a named vector
110+
# df <- dplyr::tibble(
111+
# est_df = bbmle::coef(mle_fit)[1],
112+
# est_ncp = bbmle::coef(mle_fit)[2]
113+
# )
114+
# return(df)
115+
# }
116+
#
117+
# safe_estimates <- {
118+
# purrr::possibly(
119+
# estimate_chisq_params,
120+
# otherwise = NA_real_,
121+
# quiet = TRUE
122+
# )
123+
# }
124+
#
125+
# estimates <- safe_estimates(x_term)
126+
# Define negative log-likelihood function
127+
neg_log_likelihood <- function(params) {
128+
df <- params[1]
129+
ncp <- params[2]
130+
sum_densities <- sum(dchisq(x_term, df = df, ncp = ncp, log = TRUE))
131+
return(-sum_densities)
115132
}
116133

117-
safe_estimates <- {
118-
purrr::possibly(
119-
estimate_chisq_params,
120-
otherwise = NA_real_,
121-
quiet = TRUE
122-
)
123-
}
134+
# Initial guess for parameters
135+
initial_params <- c(trunc(var(x_term)/2), trunc(mean(x_term)))
136+
137+
# Optimize parameters using optim() function
138+
opt_result <- optim(par = initial_params, fn = neg_log_likelihood)
124139

125-
estimates <- safe_estimates(x_term)
126-
doff <- estimates$est_df |> unname()
127-
ncp <- estimates$est_ncp |> unname()
140+
# Extract estimated parameters
141+
doff <- opt_result$par[1]
142+
ncp <- opt_result$par[2]
128143

129144
# Return Tibble ----
130145
if (.auto_gen_empirical) {
@@ -139,7 +154,7 @@ util_chisquare_param_estimate <- function(.x, .auto_gen_empirical = TRUE) {
139154
min = minx,
140155
max = maxx,
141156
mean = mean(x_term),
142-
degrees_of_freedom = doff,
157+
dof = doff,
143158
ncp = ncp
144159
)
145160

0 commit comments

Comments
 (0)