# R/loop.optim.R In HPbayes: Heligman Pollard mortality model parameter estimation using Bayesian Melding with Incremental Mixture Importance Sampling

#### Documented in loop.optim

```loop.optim <- function(prior, nrisk, ndeath, d = 10, theta.dim = 8, age = c(1e-05, 1, seq(5, 100, 5)))
{
lx <- nrisk
dx <- ndeath
H.k <- prior
pllwts <- prior.likewts(prior = prior, nrisk = lx, ndeath = dx)
log.like.0 <- pllwts\$log.like.0
wts.0 <- pllwts\$wts.0
B0 <- 1000 * theta.dim
q0 <- H.k
d.keep <- 0
theta.new <- H.k[wts.0 == max(wts.0), ]
keep <- H.k
ll.keep <- log.like.0
opt.mu.d <- matrix(NA, nrow = d, ncol = theta.dim)
opt.cov.d <- array(NA, dim = c(theta.dim, theta.dim, d))
prior.cov <- cov(q0)
opt.low <- apply(q0, 2, min)
opt.hi <- apply(q0, 2, max)
imp.keep <- theta.dim * 100
max.log.like.0 <- max(log.like.0)
mp8.mle <- function(theta, x.fit = age) {
p.hat <- mod8p(theta = theta, x = x.fit)
ll <- ll.binom(x = dx, n = lx, p = p.hat)
return(ll)
}
for (i in 1:d) {
out <- optim(par = theta.new, fn = mp8.mle, method = "L-BFGS-B",
lower = opt.low, upper = opt.hi, control = list(fnscale = -1,
maxit = 1e+05))
out.mu <- out\$par
if (out\$value > max.log.like.0) {
d.keep <- d.keep + 1
opt.mu.d[i, ] <- out.mu
out.hess <- hessian(func = mp8.mle, x = out\$par)
if (is.positive.definite(-out.hess)) {
out.cov <- try(solve(-out.hess))
opt.cov.d[, , i] <- out.cov
}
if (!is.positive.definite(-out.hess)) {
out.prec <- try(solve(prior.cov)) + A
if (!is.positive.definite(out.prec)) {
out.prec <- solve(prior.cov)
}
out.cov <- try(solve(out.prec))
opt.cov.d[, , i] <- out.cov
}
}

if(i==1 & out\$value <= max.log.like.0){
out.hess <- hessian(func = mp8.mle, x = out\$par)
if (is.positive.definite(-out.hess)) {
out.cov <- solve(-out.hess)
}
if (!is.positive.definite(-out.hess)) {
out.prec <- solve(prior.cov) + A
if (!is.positive.definite(out.prec)) {
out.prec <- solve(prior.cov)
}
out.cov <- solve(out.prec)
}
warning("likelihood of first local maximum does not exceed maximum 			likelihood from the prior")
}

if (i < d) {
keep <- keep[ll.keep != max(ll.keep), ]
ll.keep <- ll.keep[ll.keep != max(ll.keep)]
dist.to.mu <- mahalanobis(x = keep, center = out.mu,
cov = out.cov)
keep <- keep[rank(1/dist.to.mu) <= (d - i) * B0/d,
]
ll.keep <- ll.keep[rank(1/dist.to.mu) <= (d - i) *
B0/d]
theta.new <- keep[ll.keep == max(ll.keep), ]
}
}
return(list(opt.mu.d = opt.mu.d, opt.cov.d = opt.cov.d, theta.new = theta.new,
d.keep = d.keep, log.like.0 = log.like.0, wts.0 = wts.0))
}
```

## Try the HPbayes package in your browser

Any scripts or data that you put into this service are public.

HPbayes documentation built on May 2, 2019, 5:53 a.m.