Description Usage Arguments Value See Also Examples
The turbo
class represents results from parameter estimation in fixed-point mapping problems. The turboem
function outputs objects of class turbo
.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 | ## S3 method for class 'turbo'
print(x, ...)
## S3 method for class 'turbo'
pars(x, ...)
## S3 method for class 'turbo'
error(x, ...)
## S3 method for class 'turbo'
plot(x, which.methods = seq_along(x$method),
method.names = x$method[which.methods], xlim, ylim, ...)
## S3 method for class 'turbo'
grad(x, objfn=x$objfn, which.methods = seq_along(x$method),
method.names = x$method[which.methods], ...)
## S3 method for class 'turbo'
hessian(x, objfn=x$objfn, which.methods = seq_along(x$method),
method.names = x$method[which.methods], ...)
## S3 method for class 'turbo'
stderror(x, objfn=x$objfn, which.methods = seq_along(x$method),
method.names = x$method[which.methods], ...)
|
x |
An object of class |
which.methods |
A vector identifying for which subset of algorithms results are desired. |
method.names |
A vector of unique identifiers for the algorithms for which results are being provided. |
xlim |
Optional range for the x-axis of the trace plot. |
ylim |
Optional range for the y-axis of the trace plot. |
objfn |
Objective function. Usually this is taken to be the appropriate component of a |
... |
Additional arguments. |
|
Shows a brief summary of the results from fitting the acceleration schemes. |
|
Prints the fixed-point values across acceleration schemes at termination of the algorithms. |
|
Prints any error messages from running the acceleration schemes |
|
Shows a trace plot of the objective function value over iterations. This method is only available if the call to |
|
Calculates the gradient of the objective function evaluated at the fixed-point across acceleration schemes. Uses numerical methods from the package |
|
Calculates the Hessian of the objective function evaluated at the fixed-point across acceleration schemes. Uses numerical methods from the package |
|
Provides estimates of the standard error of the fixed-point across acceleration schemes. |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 | ###########################################################################
# Also see the vignette by typing:
# vignette("turboEM")
#
# EM algorithm for Poisson mixture estimation
fixptfn <- function(p,y) {
# The fixed point mapping giving a single E and M step of the EM algorithm
#
pnew <- rep(NA,3)
i <- 0:(length(y)-1)
zi <- p[1]*exp(-p[2])*p[2]^i / (p[1]*exp(-p[2])*p[2]^i + (1 - p[1])*exp(-p[3])*p[3]^i)
pnew[1] <- sum(y*zi)/sum(y)
pnew[2] <- sum(y*i*zi)/sum(y*zi)
pnew[3] <- sum(y*i*(1-zi))/sum(y*(1-zi))
p <- pnew
return(pnew)
}
objfn <- function(p,y) {
# Objective function whose local minimum is a fixed point
# negative log-likelihood of binary poisson mixture
i <- 0:(length(y)-1)
loglik <- y*log(p[1]*exp(-p[2])*p[2]^i/exp(lgamma(i+1)) +
(1 - p[1])*exp(-p[3])*p[3]^i/exp(lgamma(i+1)))
return ( -sum(loglik) )
}
# Real data from Hasselblad (JASA 1969)
poissmix.dat <- data.frame(death=0:9, freq=c(162,267,271,185,111,61,27,8,3,1))
y <- poissmix.dat$freq
# Use a preset seed so the example is reproducable.
require("setRNG")
old.seed <- setRNG(list(kind="Mersenne-Twister", normal.kind="Inversion",
seed=1))
p0 <- c(runif(1),runif(2,0,4)) # random starting value
# Basic EM algorithm, SQUAREM, and parabolic EM, with default settings
res1 <- turboem(par=p0, y=y, fixptfn=fixptfn, objfn=objfn, method=c("EM", "squarem", "pem"))
# Apply methods for class "turbo"
res1
pars(res1)
grad(res1)
hessian(res1)
stderror(res1)
error(res1)
# We get an error for Dynamic ECME (decme) if we do not specify the boundary function
res2 <- turboem(par=p0, y=y, fixptfn=fixptfn, objfn=objfn,
method=c("EM", "squarem", "pem", "decme"))
res2
error(res2)
# we can't plot the results, because we did not store the objective function value at each iteration
# Changing the options to store the objective function values, we can:
res1keep <- turboem(par=p0, y=y, fixptfn=fixptfn, objfn=objfn, method=c("EM", "squarem", "pem"),
control.run=list(keep.objfval=TRUE))
plot(res1keep, xlim=c(0.001, 0.02))
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.