Revision 63c8e8a453ea587001e2438a8ce51cf0e1e1675c authored by Charles J. Geyer on 23 March 2009, 00:00:00 UTC, committed by Gabor Csardi on 23 March 2009, 00:00:00 UTC
1 parent b524c08
trans.R
library(aster)
set.seed(42)
pred <- c(0, 1, 1, 2)
fam <- c(1, 1, 3, 2)
famlist <- fam.default()
nnode <- length(pred)
nind <- 7
x <- rbind(c(1, 1, 2, 4),
c(1, 1, 3, 0),
c(1, 1, 5, 2),
c(1, 0, 1, 0),
c(1, 0, 4, 0),
c(0, 0, 0, 0),
c(0, 0, 0, 0))
nrow(x) == nind
ncol(x) == nnode
root <- matrix(1, nind, nnode)
theta <- rnorm(nind * nnode, 0, .33)
theta <- matrix(theta, nind, nnode)
storage.mode(theta) <- "double"
aster:::setfam(fam.default())
out <- .C("aster_theta2phi",
nind = as.integer(nind),
nnode = as.integer(nnode),
pred = as.integer(pred),
fam = as.integer(fam),
theta = theta,
phi = matrix(as.double(0), nind, nnode))
my.phi <- theta
# poisson cumulant function
my.phi[ , 2] <- my.phi[ , 2] - exp(theta[ , 4])
# non-zero poisson cumulant function
my.phi[ , 1] <- my.phi[ , 1] - log(exp(exp(theta[ , 3])) - 1)
# bernoulli cumulant function
my.phi[ , 1] <- my.phi[ , 1] - log(1 + exp(theta[ , 2]))
all.equal(out$phi, my.phi)
tout <- .C("aster_phi2theta",
nind = as.integer(nind),
nnode = as.integer(nnode),
pred = as.integer(pred),
fam = as.integer(fam),
phi = out$phi,
theta = matrix(as.double(0), nind, nnode))
all.equal(tout$theta, theta)
storage.mode(x) <- "double"
storage.mode(root) <- "double"
cout <- .C("aster_theta2ctau",
nind = as.integer(nind),
nnode = as.integer(nnode),
pred = as.integer(pred),
fam = as.integer(fam),
theta = theta,
ctau = matrix(as.double(0), nind, nnode))
my.ctau <- theta
# poisson mean function
my.ctau[ , 4] <- exp(theta[ , 4])
# non-zero poisson mean function
foo <- exp(theta[ , 3])
my.ctau[ , 3] <- foo / (1 - exp(- foo))
# bernoulli cumulant function
my.ctau[ , 2] <- 1 / (1 + exp(- theta[ , 2]))
# bernoulli cumulant function
my.ctau[ , 1] <- 1 * 1 / (1 + exp(- theta[ , 1]))
all.equal(cout$ctau, my.ctau)
xout <- .C("aster_xpred",
nind = as.integer(nind),
nnode = as.integer(nnode),
pred = as.integer(pred),
fam = as.integer(fam),
x = x,
root = root,
xpred = matrix(as.double(0), nind, nnode))
my.xpred <- x
my.xpred[ , 1] <- root[ , 1]
my.xpred[ , 2] <- x[ , 1]
my.xpred[ , 3] <- x[ , 1]
my.xpred[ , 4] <- x[ , 2]
all.equal(xout$xpred, my.xpred)
tout <- .C("aster_ctau2tau",
nind = as.integer(nind),
nnode = as.integer(nnode),
pred = as.integer(pred),
fam = as.integer(fam),
root = root,
ctau = cout$ctau,
tau = matrix(as.double(0), nind, nnode))
my.tau <- my.ctau
my.tau[ , 1] <- my.xpred[ , 1] * my.ctau[ , 1]
my.tau[ , 2] <- my.tau[ , 1] * my.ctau[ , 2]
my.tau[ , 3] <- my.tau[ , 1] * my.ctau[ , 3]
my.tau[ , 4] <- my.tau[ , 2] * my.ctau[ , 4]
all.equal(tout$tau, my.tau)
##########
wout <- .C("aster_theta2whatsis",
nind = as.integer(nind),
nnode = as.integer(nnode),
pred = as.integer(pred),
fam = as.integer(fam),
deriv = as.integer(0),
theta = as.double(theta),
result = matrix(as.double(0), nind, nnode))
my.what <- NaN * wout$result
for (i in 1:nind)
for (j in 1:nnode) {
ifam <- fam[j]
my.what[i, j] <- famfun(famlist[[ifam]], 0, theta[i, j])
}
all.equal(wout$result, my.what)
##########
aster:::setfam(fam.default())
wout <- .C("aster_theta2whatsis",
nind = as.integer(nind),
nnode = as.integer(nnode),
pred = as.integer(pred),
fam = as.integer(fam),
deriv = as.integer(1),
theta = as.double(theta),
result = matrix(as.double(0), nind, nnode))
my.what <- NaN * wout$result
for (i in 1:nind)
for (j in 1:nnode) {
ifam <- fam[j]
my.what[i, j] <- famfun(famlist[[ifam]], 1, theta[i, j])
}
all.equal(wout$result, my.what)
##########
aster:::setfam(fam.default())
wout <- .C("aster_theta2whatsis",
nind = as.integer(nind),
nnode = as.integer(nnode),
pred = as.integer(pred),
fam = as.integer(fam),
deriv = as.integer(2),
theta = as.double(theta),
result = matrix(as.double(0), nind, nnode))
my.what <- NaN * wout$result
for (i in 1:nind)
for (j in 1:nnode) {
ifam <- fam[j]
my.what[i, j] <- famfun(famlist[[ifam]], 2, theta[i, j])
}
all.equal(wout$result, my.what)
##########
aster:::setfam(fam.default())
vout <- .C("aster_tt2var",
nind = as.integer(nind),
nnode = as.integer(nnode),
pred = as.integer(pred),
fam = as.integer(fam),
x = as.double(x),
root = as.double(root),
theta = as.double(theta),
tau = as.double(tout$tau),
var = matrix(as.double(0), nind, nnode))
my.var <- NaN * vout$var
for (i in 1:nind)
for (j in 1:nnode) {
ifam <- fam[j]
thefam <- famlist[[ifam]]
k <- pred[j]
if (k > 0) {
my.var[i, j] <- famfun(thefam, 2, theta[i, j]) * tout$tau[i, k] +
famfun(thefam, 1, theta[i, j])^2 * my.var[i, k]
} else {
my.var[i, j] <- famfun(thefam, 2, theta[i, j]) * xout$xpred[i, j]
}
}
all.equal(vout$var, my.var)
Computing file changes ...