Nothing
library(trust)
options(digits = 3)
##### four-way contingency table with all two-way interactions
d <- c(3, 4, 5, 6)
n <- 1000
##### model matrix
m <- NULL
for (i in 1:d[1]) {
for (j in 1:d[2]) {
mfoo <- array(0, dim = d)
mfoo[i, j, , ] <- 1
mfoo <- as.vector(mfoo)
m <- cbind(m, mfoo)
}
}
for (i in 1:d[1]) {
for (j in 1:d[3]) {
mfoo <- array(0, dim = d)
mfoo[i, , j, ] <- 1
mfoo <- as.vector(mfoo)
m <- cbind(m, mfoo)
}
}
for (i in 1:d[1]) {
for (j in 1:d[4]) {
mfoo <- array(0, dim = d)
mfoo[i, , , j] <- 1
mfoo <- as.vector(mfoo)
m <- cbind(m, mfoo)
}
}
for (i in 1:d[2]) {
for (j in 1:d[3]) {
mfoo <- array(0, dim = d)
mfoo[ , i, j, ] <- 1
mfoo <- as.vector(mfoo)
m <- cbind(m, mfoo)
}
}
for (i in 1:d[2]) {
for (j in 1:d[4]) {
mfoo <- array(0, dim = d)
mfoo[ , i, , j] <- 1
mfoo <- as.vector(mfoo)
m <- cbind(m, mfoo)
}
}
for (i in 1:d[3]) {
for (j in 1:d[4]) {
mfoo <- array(0, dim = d)
mfoo[ , , i, j] <- 1
mfoo <- as.vector(mfoo)
m <- cbind(m, mfoo)
}
}
dimnames(m) <- NULL
foo <- qr(m)
m <- m[ , foo$pivot[seq(1, foo$rank)]]
##### true parameter value
set.seed(42)
theta.true <- 0.25 * rnorm(ncol(m))
theta.true <- round(theta.true, 5)
##### simulate data
eta <- as.numeric(m %*% theta.true)
p <- exp(eta)
p <- p / sum(p)
x <- sample(nrow(m), n, replace = TRUE, prob = p)
x <- tabulate(x, nbins = nrow(m))
##### save data
iffy <- try(read.table("fred.txt"), silent = TRUE)
if (inherits(iffy, "try-error")) {
data <- data.frame(x = x, m = m)
write.table(data, file = "fred.txt", row.names = FALSE)
}
data <- read.table(file = "fred.txt", header = TRUE)
x <- data$x
data$x <- NULL
m <- as.matrix(data)
dimnames(m) <- NULL
##### log likelihood
objfun <- function(theta) {
eta <- as.numeric(m %*% theta)
p <- exp(eta)
f <- sum(x * eta - p)
g <- as.numeric(t(x - p) %*% m)
B <- sweep(- m, 1, p, "*")
B <- t(m) %*% B
list(value = f, gradient = g, hessian = B)
}
##### check it
sally <- objfun(theta.true)
epsilon <- 1e-8
mygrad <- double(length(theta.true))
for (i in 1:length(mygrad)) {
theta.eps <- theta.true
theta.eps[i] <- theta.true[i] + epsilon
sally.eps <- objfun(theta.eps)
mygrad[i] <- (sally.eps$value - sally$value) / epsilon
}
all.equal(sally$gradient, mygrad, tolerance = length(mygrad) * epsilon)
myhess <- matrix(NA, length(theta.true), length(theta.true))
for (i in 1:length(mygrad)) {
theta.eps <- theta.true
theta.eps[i] <- theta.true[i] + epsilon
sally.eps <- objfun(theta.eps)
myhess[i, ] <- (sally.eps$gradient - sally$gradient) / epsilon
}
all.equal(sally$hessian, myhess, tolerance = length(mygrad) * epsilon)
fred <- trust(objfun, theta.true, 1, sqrt(ncol(m)), minimize = FALSE)
fran <- glm.fit(m, x, family = poisson(), intercept = FALSE)
all.equal(fran$coefficients, fred$argument)
fred <- trust(objfun, rep(0, length(theta.true)), 1, sqrt(ncol(m)),
minimize = FALSE, blather = TRUE)
fred$converged
#### CRAN don't like: gives different results on different hardware
#### or different compiler flags
## ceiling(log10(max(abs(fred$gradient))))
## length(fred$r)
## data.frame(type = fred$steptype, rho = fred$rho, change = fred$preddiff,
## accept = fred$accept, r = fred$r)
## (fred$stepnorm / fred$r)[fred$accept & fred$steptype != "Newton"]
fred <- trust(objfun, rep(-5, length(theta.true)), 1, sqrt(ncol(m)),
minimize = FALSE, blather = TRUE)
fred$converged
#### CRAN don't like: gives different results on different hardware
#### or different compiler flags
## ceiling(log10(max(abs(fred$gradient))))
## length(fred$r)
## data.frame(type = fred$steptype, rho = fred$rho, change = fred$preddiff,
## accept = fred$accept, r = fred$r)
## (fred$stepnorm / fred$r)[fred$accept & fred$steptype != "Newton"]
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.