#
#mom_and_ll <- function(w, Ve, k, p_n){
# mom_est <- mom_ma_gamma(w, Ve, mean ( k ) * (1-p_n) )
# shape_mom <- mom_est[["shape"]]
# rate_mom <- mom_est[["rate"]]
#
# LL <- dma_gamma_known(w, Ve=Ve, shape=shape_mom, rate=rate_mom, p_neutral = p_n, k=k)
# c( LL_mom=LL, p_n=p_n, shape_mom = shape_mom, rate_mom = rate_mom )
#}
#
#
#mom_and_maxamize <- function(w, Ve, k, p_n){
# res <- mom_and_ll(w, Ve, k, p_n)
# mod <- fit_gamma_known(w,
# start=list(shape=res[["shape_mom"]], rate=res[["rate_mom"]]),
# fixed=list(Ve=Ve, p_neutral = p_n), k)
# c(res, mod@coef, maxLL=-mod@min)
#
#}
#
#
#AIC <- function(ests) 6- 2 *ests[5,]
#ests <- sapply(Ut, mom_and_ll, w=w, Ve=0.001)
#Are all of the requred args either set or give a starting value
check_args <- function(required){
parent_args <- sys.frame(sys.parent())
all_args <- (c(names(parent_args$start), names(parent_args$fixed)))
if( !all (required %in% all_args )) {
missing <- paste(required[!(required %in% all_args)])
msg <- paste("Must set fixed or starting values for the following parameters:\n",
paste(missing, collapse=" "))
stop(msg, call.=FALSE)
}
}
#Create contstraints for parameters that need them
process_constraints <- function(starting_vals, constraint, lower=TRUE){
bounds <- constraint[names(starting_vals)]
NO_CONSTRAINT <- if(lower) -Inf else Inf
res <- vapply(bounds, function(x) if(is.null(x)) NO_CONSTRAINT else x, FUN.VALUE=0.0)
if(all(is.infinite(res))){
return(NULL)
}
unname(res)
}
# closures to generating fitting functions
make_dfe_fitting_fxn <- function(Q, all_params, verbose=TRUE, lower=list(), upper=list()){
f <- function(obs, start=list(), fixed=list(), verbose=versbose, ...){
check_args(all_params)
#Set up the specific fitting function (with fixed values given and (to
#match optimx) variables given as a single vector)
Q_args <- formals(Q)
Q_args$obs <- obs
Q_args[names(fixed)] <- fixed
Q_args[names(start)] <- NULL
#Are we using contraints? If so, let's set them up
lower_bound <- process_constraints(start, lower, lower=TRUE)
upper_bound <- process_constraints(start, upper, lower=FALSE)
LL <- function(x) {
res <- -do.call(Q, c(x, Q_args))
if (is.nan(res)){
return( .Machine$double.xmax )
}
if (is.infite(res)){
.dfe_bad_obs <<- x
.dfe_bad_args <<- Q_args
return (.Machine$double.xmax )
}
res
}
if(is.null(upper_bound) & is.null(lower_bound)){
return( optimx(unlist(start), LL, ...) )
}
LL_args <- c(list(method="L-BFGS-B", par=unlist(start), fn=LL), list(...))
if(!is.null(lower_bound)){
LL_args$lower <- lower_bound
}
if(!is.null(upper_bound)){
LL_args$upper <- upper_bound
}
do.call(optimx, LL_args)
}
f
}
# @export
#fit_gamma_known <- make_dfe_fitting_fxn(dma_gamma_known,
# c("shape", "rate", "Ve", "p_neutral", "k"),
# lower=list(shape=1e-8, rate=1e-8, Ve=1e-8, p_neutral=0),
# upper = list(p_neutral=1.0))
#mom_lik_sim <- function(shape, rate, Ut, pn, n, Ve=0.01){
# pn_grid <- seq(0, 1, 0.025)
# mu <- rpois(n,Ut)
# w <- rma_known_gamma(shape=shape, rate=rate, Ve=Ve, k=mu, p_neutral=pn)
# res <- vapply(pn_grid, mom_and_ll, w=w, Ve=0.01, k=mu, FUN.VALUE=rep(0.01,4))
# simmed_pn <- which(pn_grid == pn)
# winner <- which.max(res[1,])
# delta_LL <- unname(res[1,winner] - res[1,simmed_pn])
# params = as.list(match.call())[-1]
# unlist(c(params, LL_dif=delta_LL, res[,winner]))
#
#}
#
#mod <- fit_gamma_known(w,
# fixed=list(Ve=0.01, shape=1, rate=10, k=mu),
# start=list(p_neutral = 0.1) )
#
#
#
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.