stats_function(eq_c = NULL, eq_d = NULL, par_c = NULL, par_d = NULL, fixed = 0)

eq_d <- c("ASC1 * 1 + B11_dur * dur_1" , "ASC2 * 1 + B12_dur * dur_2",
"ASC3 * 1 + B13_dur * dur_3")
eq_c <- c("Tw ~ tw*w + ph1*Tc", "Tf1 ~ (1+w)^tw + ph1^3*Tc")
parl <- c("tw", "ph1")
par_d <- c(paste0("ASC", 1:3), paste0("B1", 1:3, "_dur"))
stfunc <- stats_function(eq_c, eq_d, parl,par_d, fixed=3)
data <- matrix(runif(1000, min=0.001, max=50), ncol=8)
data <- data.frame(data)
names(data) <- c("dur_1", "dur_2", "dur_3", "w", "Tc", "avl_1", "avl_2", "avl_3")
parv <- c(0.5, 1, 1.5, 2, 1, -0.3, 0.2, -0.8)
methodopt <- "BHHH"
separatenmm <- TRUE
env <- environment()
fnames <- c("prob_func", "cont_func")
sapply(1:(length(stfunc)-1), function(x)assign(fnames[x], stfunc[[x]], envir=env))
eval(parse(text=paste0("environment(", fnames[1:(length(stfunc)-1)], ") <- env")))
probs <- prob_func(parv)
apply(probs, 2, mean)
cont <- cont_func(parv)
apply(cont, 2, mean)