- dev <- function(a, r)
- {
- -2 * (lfactorial(N - 1)
- + (N - 2) * log(r)
- + (3 * N) * log(1 - a)
- + 2 * r * sum(x[2:N])
- - 2 * sum(log(exp(r * x[2:N]) - a))
- + r * sum(trm.br)
- + sum((S - 1) * log(exp(r * trm.br) - 1))
- - sum((S + 1) * log(exp(r * trm.br) - a)))
- }
- out <- nlm(function(p) dev(p[1], p[2]), c(0, 0.2), hessian = TRUE)
- if (out$estimate[1] < 0) {
- out <- nlm(function(p) dev(0, p), 0.2, hessian = TRUE)
- para <- c(0, out$estimate)
- se <- c(0, sqrt(diag(solve(out$hessian))))
- }
- else {
- para <- out$estimate
- se <- sqrt(diag(solve(out$hessian)))
+ if (conditional) {
+ dev <- function(a, r) {
+ if (a >= 1 || a < 0 || r <= 0) return(1e50)
+ ert <- exp(r * trm.br)
+ zeta <- (ert - 1)/(ert - a)
+ -2 * (lfactorial(N - 1)
+ + (N - 2) * log(r)
+ + N * log(1 - a)
+ + 2 * r * sum(x[2:N])
+ - 2 * sum(log(exp(r * x[2:N]) - a))
+ + sum(log(1 - zeta) + (S - 1)*log(zeta)))
+ }
+ } else {
+ dev <- function(a, r) {
+ if (a >= 1 || a < 0 || r <= 0) return(1e50)
+ -2 * (lfactorial(N - 1)
+ + (N - 2) * log(r)
+ + (3 * N) * log(1 - a)
+ + 2 * r * sum(x[2:N])
+ - 2 * sum(log(exp(r * x[2:N]) - a))
+ + r * sum(trm.br)
+ + sum((S - 1) * log(exp(r * trm.br) - 1))
+ - sum((S + 1) * log(exp(r * trm.br) - a)))
+ }