Nothing
library(aster2)
set.seed(42)
data(test1)
fred <- asterdata(test1,
vars = c("m1", "m2", "m3", "n1", "n2", "b1", "p1", "z1"),
pred = c(0, 0, 0, 1, 1, 2, 3, 6), group = c(0, 1, 2, 0, 4, 0, 0, 0),
code = c(1, 1, 1, 2, 2, 3, 4, 5),
families = list(fam.multinomial(3), "normal.location.scale",
"bernoulli", "poisson", "zero.truncated.poisson"))
########## valid theta ##########
theta <- rnorm(nrow(fred$redata)) * 0.1
try(validtheta(fred, theta))
must.be.negative <- as.character(fred$redata$varb) == "n2"
theta[must.be.negative] <- (- abs(theta[must.be.negative]))
is.validtheta(fred, theta)
########## theta to phi ##########
phi <- transformSaturated(theta, fred, from = "theta", to = "phi")
group <- fred$regroup
group.idx <- seq(along = group)
group.grp <- cbind(group.idx)
repeat {
group.idx <- c(0, group)[group.idx + 1]
if (all(group.idx == 0)) break
group.grp <- cbind(group.grp, group.idx)
}
outies <- sort(unique(as.vector(group.grp[ , -1])))
outies <- outies[outies != 0]
if (length(outies) > 0)
group.grp <- group.grp[- outies, ]
group.grp <- as.data.frame(t(group.grp))
group.grp <- as.list(group.grp)
names(group.grp) <- NULL
group.grp <- lapply(group.grp, function(x) x[x != 0])
group.grp <- lapply(group.grp, rev)
identical(as.numeric(seq(along = group)), as.numeric(sort(unlist(group.grp))))
delta <- fred$redelta
families <- fred$families[fred$recode]
z <- sapply(group.grp, function(i) cumulant(theta[i],
families[[i[1]]], delta = delta[i])$zeroth)
group.pred <- fred$repred[sapply(group.grp, function(x) x[1])]
myphi <- theta
for (i in seq(along = group.grp)) {
j <- group.pred[i]
if (j != 0) {
myphi[j] <- myphi[j] - z[i]
}
}
all.equal(phi, myphi)
########## theta to phi derivative ##########
dtheta <- rnorm(length(theta))
dphi <- transformSaturated(theta, fred, from = "theta", to = "phi",
differential = dtheta)
xi <- lapply(group.grp, function(i) cumulant(theta[i],
families[[i[1]]], delta = delta[i], deriv = 1)$first)
mydphi <- dtheta
for (i in rev(seq(along = group.grp))) {
j <- group.pred[i]
if (j != 0) {
mydphi[j] <- mydphi[j] - sum(xi[[i]] * dtheta[group.grp[[i]]])
}
}
all.equal(dphi, mydphi)
epsilon <- 1e-8
mymydphi <- (transformSaturated(theta + epsilon * dtheta, fred,
from = "theta", to = "phi") -
transformSaturated(theta - epsilon * dtheta, fred,
from = "theta", to = "phi")) / (2 * epsilon)
all.equal(dphi, mymydphi)
########## phi to theta ##########
newtheta <- transformSaturated(phi, fred, from = "phi", to = "theta")
all.equal(theta, newtheta)
########## dphi to dtheta ##########
newdtheta <- transformSaturated(phi, fred, from = "phi", to = "theta",
differential = dphi)
all.equal(dtheta, newdtheta)
########## theta to xi ##########
xi <- transformSaturated(theta, fred, from = "theta", to = "xi")
is.validxi(fred, xi)
myxi <- rep(NA, length(xi))
for (i in seq(along = group.grp)) {
j = group.grp[[i]]
out <- cumulant(theta[j], families[[j[1]]], delta = delta[j], deriv = 1)
myxi[j] <- out$first
}
all.equal(xi, myxi)
########## dtheta to dxi ##########
dxi <- transformSaturated(theta, fred, from = "theta", to = "xi",
differential = dtheta)
mydxi <- rep(NA, length(dxi))
for (i in seq(along = group.grp)) {
j = group.grp[[i]]
out <- cumulant(theta[j], families[[j[1]]], delta = delta[j], deriv = 2)
mydxi[j] <- as.numeric(out$second %*% dtheta[j])
}
all.equal(dxi, mydxi)
########## xi to theta ##########
mytheta <- transformSaturated(xi, fred, from = "xi", to = "theta")
is.same(theta, mytheta, fred, parm.type = "theta")
## mymyxi <- transformSaturated(mytheta, fred, from = "theta", to = "xi")
## all.equal(xi, mymyxi)
########## dxi to dtheta ##########
mydtheta <- transformSaturated(xi, fred, from = "xi", to = "theta",
differential = dxi)
mymydtheta <- (
transformSaturated(xi + epsilon * dxi, fred, from = "xi", to = "theta") -
transformSaturated(xi - epsilon * dxi, fred, from = "xi", to = "theta")
) / (2 * epsilon)
all.equal(mydtheta, mymydtheta, tol = 1e-7)
is.same(dtheta, mydtheta, fred, parm.type = "theta")
########## xi to mu ##########
mu <- transformSaturated(xi, fred, from = "xi", to = "mu")
mymu <- rep(NA, length(mu))
while (any(is.na(mymu))) {
mymu.pred <- ifelse(fred$repred == 0, fred$initial,
c(0, mymu)[fred$repred + 1])
mymu <- xi * mymu.pred
}
all.equal(mu, mymu)
## mupred <- ifelse(fred$repred == 0, fred$initial, c(0, mu)[fred$repred + 1])
## mymymyxi <- mu / mupred
## all.equal(xi, mymymyxi)
########## dxi to dmu ##########
dmu <- transformSaturated(xi, fred, from = "xi", to = "mu",
differential = dxi)
mydmu <- (
transformSaturated(xi + epsilon * dxi, fred, from = "xi", to = "mu") -
transformSaturated(xi - epsilon * dxi, fred, from = "xi", to = "mu")
) / (2 * epsilon)
all.equal(dmu, mydmu)
########## mu to xi ##########
mymymymyxi <- transformSaturated(mu, fred, from = "mu", to = "xi")
all.equal(xi, mymymymyxi)
########## dmu to dxi ##########
mymymymydxi <- transformSaturated(mu, fred, from = "mu", to = "xi",
differential = dmu)
all.equal(dxi, mymymymydxi)
########## theta to mu ##########
mu.foo <- transformSaturated(theta, fred, from = "theta", to = "mu")
all.equal(mu, mu.foo)
########## dtheta to dmu ##########
dmu.foo <- transformSaturated(theta, fred, from = "theta", to = "mu",
differential = dtheta)
all.equal(dmu, dmu.foo)
########## xi to phi ##########
phi.foo <- transformSaturated(xi, fred, from = "xi", to = "phi")
is.same(phi, phi.foo, fred, parm.type ="phi")
########## dxi to dphi ##########
dphi.foo <- transformSaturated(xi, fred, from = "xi", to = "phi",
differential = dxi)
is.same(dphi, dphi.foo, fred, parm.type ="phi")
########## phi to xi ##########
xi.foo <- transformSaturated(phi, fred, from = "phi", to = "xi")
all.equal(xi, xi.foo)
########## dphi to dxi ##########
dxi.foo <- transformSaturated(phi, fred, from = "phi", to = "xi",
differential = dphi)
all.equal(dxi, dxi.foo)
########## phi to mu ##########
mu.foo <- transformSaturated(phi, fred, from = "phi", to = "mu")
all.equal(mu, mu.foo)
########## dphi to dmu ##########
dmu.foo <- transformSaturated(phi, fred, from = "phi", to = "mu",
differential = dphi)
all.equal(dmu, dmu.foo)
########## mu to theta ##########
theta.bar <- transformSaturated(mu, fred, from = "mu", to = "theta")
is.same(theta, theta.bar, fred, parm.type ="theta")
########## dmu to dtheta ##########
dtheta.bar <- transformSaturated(mu, fred, from = "mu", to = "theta",
differential = dmu)
is.same(dtheta, dtheta.bar, fred, parm.type ="theta")
########## mu to phi ##########
phi.bar <- transformSaturated(mu, fred, from = "mu", to = "phi")
is.same(phi, phi.bar, fred, parm.type ="phi")
########## dmu to dphi ##########
dphi.bar <- transformSaturated(mu, fred, from = "mu", to = "phi",
differential = dmu)
is.same(dphi, dphi.bar, fred, parm.type ="phi")
########## valid xi ##########
xi <- rnorm(nrow(fred$redata))
try(validxi(fred, xi))
nind <- nrow(fred$redata) / nlevels(fred$redata$varb)
xi[as.character(fred$redata$varb) == "z1"] <- 1 + rexp(nind)
try(validxi(fred, xi))
xi[as.character(fred$redata$varb) == "p1"] <- rexp(nind)
try(validxi(fred, xi))
xi[as.character(fred$redata$varb) == "b1"] <- runif(nind)
try(validxi(fred, xi))
xi[as.character(fred$redata$varb) == "n2"] <-
xi[as.character(fred$redata$varb) == "n1"]^2 + rchisq(nind, df = 1)
try(validxi(fred, xi))
xi[as.character(fred$redata$varb) == "m1"] <- rexp(nind)
xi[as.character(fred$redata$varb) == "m2"] <- rexp(nind)
xi[as.character(fred$redata$varb) == "m3"] <- rexp(nind)
try(validxi(fred, xi))
thesum <- xi[as.character(fred$redata$varb) == "m1"] +
xi[as.character(fred$redata$varb) == "m2"] +
xi[as.character(fred$redata$varb) == "m3"]
xi[as.character(fred$redata$varb) == "m1"] <-
xi[as.character(fred$redata$varb) == "m1"] / thesum
xi[as.character(fred$redata$varb) == "m2"] <-
xi[as.character(fred$redata$varb) == "m2"] / thesum
xi[as.character(fred$redata$varb) == "m3"] <-
xi[as.character(fred$redata$varb) == "m3"] / thesum
is.validxi(fred, xi)
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.