Nothing
### R code from vignette source 'SQUAREM.Rnw'
###################################################
### code chunk number 1: load
###################################################
library("SQUAREM")
###################################################
### code chunk number 2: help
###################################################
help(package=SQUAREM)
###################################################
### code chunk number 3: rng
###################################################
require("setRNG")
setRNG(list(kind="Wichmann-Hill", normal.kind="Box-Muller", seed=123))
###################################################
### code chunk number 4: data
###################################################
poissmix.dat <- data.frame(death=0:9, freq=c(162,267,271,185,111,61,27,8,3,1))
###################################################
### code chunk number 5: init
###################################################
y <- poissmix.dat$freq
tol <- 1.e-08
setRNG(list(kind="Wichmann-Hill", normal.kind="Box-Muller", seed=123))
p0 <- c(runif(1),runif(2,0,6))
###################################################
### code chunk number 6: poissmix
###################################################
poissmix.em <- function(p,y) {
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)
}
###################################################
### code chunk number 7: loglik
###################################################
poissmix.loglik <- function(p,y) {
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) )
}
###################################################
### code chunk number 8: em
###################################################
pf1 <- fpiter(p=p0, y=y, fixptfn=poissmix.em, objfn=poissmix.loglik,
control=list(tol=tol))
pf1
###################################################
### code chunk number 9: squarem
###################################################
pf2 <- squarem(p=p0, y=y, fixptfn=poissmix.em, objfn=poissmix.loglik,
control=list(tol=tol))
pf2
###################################################
### code chunk number 10: squarem2
###################################################
pf3 <- squarem(p=p0, y=y, fixptfn=poissmix.em, control=list(tol=tol))
pf3
###################################################
### code chunk number 11: power
###################################################
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
}
###################################################
### code chunk number 12: bodewig
###################################################
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)
###################################################
### code chunk number 13: accelerate
###################################################
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))
# dominant eigenvector
ans1
# dominant eigenvalue
c(t(ans1$par) %*% bodewig.mat %*% ans1$par) / c(crossprod(ans1$par))
###################################################
### code chunk number 14: sq.bodewig
###################################################
ans2 <- squarem(p0, fixptfn=power.method, A=bodewig.mat)
ans2
ans2$par <- ans2$par / sqrt(sum(ans2$par^2))
c(t(ans2$par) %*% bodewig.mat %*% ans2$par) / c(crossprod(ans2$par))
###################################################
### code chunk number 15: sq2.bodewig
###################################################
ans3 <- squarem(p0, fixptfn=power.method, A=bodewig.mat, control=list(step.min0 = 0.5))
ans3
ans3$par <- ans3$par / sqrt(sum(ans3$par^2))
# eigenvalue
c(t(ans3$par) %*% bodewig.mat %*% ans3$par) / c(crossprod(ans3$par))
###################################################
### code chunk number 16: sq3.bodewig
###################################################
# Third-order SQUAREM
ans4 <- squarem(p0, fixptfn=power.method, A=bodewig.mat, control=list(K=3, method="rre"))
ans4
ans4$par <- ans4$par / sqrt(sum(ans4$par^2))
# eigenvalue
c(t(ans4$par) %*% bodewig.mat %*% ans4$par) / c(crossprod(ans4$par))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.