fpiter.RdA function to implement the fixed-point iteration algorithm. This includes monotone, contraction mappings including EM and MM algorithms
fpiter(par, fixptfn, objfn=NULL, control=list( ), ...)A vector of parameters denoting the initial guess for the fixed-point iteration.
A vector function, \(F\) that denotes the fixed-point mapping. This function is the most essential input in the package. It should accept a parameter vector as input and should return a parameter vector of same length. This function defines the fixed-point iteration: \(x_{k+1} = F(x_k)\). In the case of EM algorithm, \(F\) defines a single E and M step.
This is a scalar function, $L$, that denotes a ”merit”
function which attains its local minimum at the fixed-point of $F$.
This function should accept a parameter vector as input and should
return a scalar value. In the EM algorithm, the merit function \(L\)
is the log-likelihood. In some problems, a natural merit function may
not exist, in which case the algorithm works with only fixptfn.
The merit function function objfn does not have to be specified,
even when a natural merit function is available, especially when its
computation is expensive.
A list of control parameters to pass on to the algorithm. Full names of control list elements must be specified, otherwise, user-specifications are ignored. See *Details* below.
Arguments passed to fixptfn and objfn.
A list with the following components:
Parameter,\(x*\) that are the fixed-point of \(F\) such that \(x* = F(x*)\), if convergence is successful.
The value of the objective function \(L\) at termination.
Number of times the fixed-point function fixptfn was evaluated.
Number of times the objective function objfn was evaluated.
An integer code indicating type of convergence.
0 indicates successful convergence,
whereas 1 denotes failure to converge.
control is list of control parameters for the algorithm.
1.e-07.
maxiter - an integer denoting the maximum limit on the number of
evaluations of fixptfn, \(F\). Default is 1500.
trace - a logical variable denoting whether some of the intermediate
results of iterations should be displayed to the user.
Default is FALSE.
intermed
- a logical variable denoting whether the intermediate
results of iterations should be returned. If set to TRUE, the function
will return a matrix where each row corresponds to parameters at each iteration,
along with the corresponding fixed-point residual value. Default is FALSE.
R Varadhan and C Roland (2008), Simple and globally convergent numerical schemes for accelerating the convergence of any EM algorithm, Scandinavian Journal of Statistics, 35:335-353.
C Roland, R Varadhan, and CE Frangakis (2007), Squared polynomial extrapolation methods with cycling: an application to the positron emission tomography problem, Numerical Algorithms, 44:159-172.
##############################################################################
# Example 1: EM algorithm for Poisson mixture estimation
poissmix.em <- 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)
}
poissmix.loglik <- 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
tol <- 1.e-08
# Use a preset seed so the example is reproducible.
require("setRNG")
#> Loading required package: setRNG
old.seed <- setRNG(list(kind="Mersenne-Twister", normal.kind="Inversion",
seed=54321))
p0 <- c(runif(1),runif(2,0,4)) # random starting value
# Basic EM algorithm
pf1 <- fpiter(p=p0, y=y, fixptfn=poissmix.em, objfn=poissmix.loglik, control=list(tol=tol))
##############################################################################
# Example 2: Accelerating the convergence of power method iteration for
# finding the dominant eigenvector of a matrix
power.method <- function(x, A) {
# Defines one iteration of the power method
# x = starting guess for dominant eigenvector
# A = a square matrix
ax <- as.numeric(A %*% x)
f <- ax / sqrt(as.numeric(crossprod(ax)))
f
}
# Finding the dominant eigenvector of the Bodewig matrix
# This is a famous matrix for which power method has trouble converging
# See, for example, Sidi, Ford, and Smith (SIAM Review, 1988)
#
# Note: there are two eigenvalues that are equally dominant,
# but have opposite signs.
# Sometimes the power method finds the eigenvector corresponding to the
# large positive eigenvalue, but other times it finds the eigenvector
# corresponding to the large negative eigenvalue
b <- c(2, 1, 3, 4, 1, -3, 1, 5, 3, 1, 6, -2, 4, 5, -2, -1)
bodewig.mat <- matrix(b,4,4)
eigen(bodewig.mat)
#> eigen() decomposition
#> $values
#> [1] 7.932905 5.668864 -1.573191 -8.028578
#>
#> $vectors
#> [,1] [,2] [,3] [,4]
#> [1,] -0.5601445 0.3787027 0.6880479 0.2634624
#> [2,] -0.2116328 0.3624190 -0.6241229 0.6590407
#> [3,] -0.7767083 -0.5379352 -0.2598009 -0.1996335
#> [4,] -0.1953816 0.6601988 -0.2637503 -0.6755734
#>
p0 <- rnorm(4)
# Standard power method iteration
ans1 <- fpiter(p0, fixptfn=power.method, A=bodewig.mat)
# re-scaling the eigenvector so that it has unit length
ans1$par <- ans1$par / sqrt(sum(ans1$par^2))
ans1
#> $par
#> [1] -0.2634624 -0.6590407 0.1996335 0.6755734
#>
#> $value.objfn
#> [1] NA
#>
#> $fpevals
#> [1] 5000
#>
#> $objfevals
#> [1] 0
#>
#> $convergence
#> [1] FALSE
#>