data("sexratio")
model <- SSModel(Male ~ SSMtrend(1, Q = list(NA)), u = sexratio[,"Total"], data = sexratio,
distribution = "binomial")
fit <- fitSSM(model, inits = -15, method = "BFGS")
fit$model$Q #1.107652e-06
# Computing confidence intervals for sex ratio
# Uses importance sampling on response scale (1000 samples with antithetics)
set.seed(1)
imp <- importanceSSM(fit$model, nsim = 250, antithetics = TRUE)
sexratio.smooth <- numeric(length(model$y))
sexratio.ci <- matrix(0, length(model$y), 2)
w <- imp$w/sum(imp$w)
for(i in 1:length(model$y)){
sexr <- exp(imp$sample[i,1,])
sexratio.smooth[i]<-sum(sexr*w)
oo <- order(sexr)
sexratio.ci[i,] <- c(sexr[oo][which.min(abs(cumsum(w[oo]) - 0.05))],
sexr[oo][which.min(abs(cumsum(w[oo]) - 0.95))])
}
if (FALSE) {
# Filtered estimates
impf <- importanceSSM(fit$model, nsim = 250, antithetics = TRUE,filtered=TRUE)
sexratio.filter <- rep(NA,length(model$y))
sexratio.fci <- matrix(NA, length(model$y), 2)
w <- impf$w/rowSums(impf$w)
for(i in 2:length(model$y)){
sexr <- exp(impf$sample[i,1,])
sexratio.filter[i] <- sum(sexr*w[i,])
oo<-order(sexr)
sexratio.fci[i,] <- c(sexr[oo][which.min(abs(cumsum(w[i,oo]) - 0.05))],
sexr[oo][which.min(abs(cumsum(w[i,oo]) - 0.95))])
}
ts.plot(cbind(sexratio.smooth,sexratio.ci,sexratio.filter,sexratio.fci),
col=c(1,1,1,2,2,2),lty=c(1,2,2,1,2,2))
}
Run the code above in your browser using DataLab