Nothing
################################################################
## LINKS FOR DOUBLE METHODS:
################################################################
doubleduotrio <- function() {
doubleduotrio <- binomial()
doubleduotrio$link <- "Link for the duo-trio double test"
doubleduotrio$linkinv <- function(eta) {
ok <- eta > 0 & eta < 20
eta[eta <= 0] <- 0.25
eta[eta >= 20] <- 1
if(sum(ok)) {
eta.ok <- eta[ok]
pnorm.eta.2 <- pnorm(eta.ok * sqrt(1/2))
pnorm.eta.6 <- pnorm(eta.ok * sqrt(1/6))
eta[ok] <-
(1 - pnorm.eta.2 - pnorm.eta.6 + 2 * pnorm.eta.2 * pnorm.eta.6)^2
}
pmin(pmax(eta, 0.25), 1) ## restrict to [0.25, 1] - just to be sure
}
doubleduotrio$mu.eta <- function(eta) {
ok <- eta > 0 ## no upper limit
eta[!ok] <- 0
eta[eta >= 20] <- 1 ## should we include this restriction?
if(sum(ok)) {
eta.ok <- eta[ok]
pnorm.eta.2 <- pnorm(eta.ok * sqrt(1/2))
pnorm.eta.6 <- pnorm(eta.ok * sqrt(1/6))
sqrt.2 <- sqrt(1/2)
sqrt.6 <- sqrt(1/6)
eta.2 <- eta.ok * sqrt.2
eta.6 <- eta.ok * sqrt.6
A <- dnorm(eta.2) * sqrt.2
B <- dnorm(eta.6) * sqrt.6
C <- dnorm(eta.2)
D <- pnorm(eta.6) * sqrt.2
E <- pnorm(eta.2)
eta[ok] <- 2 * (1 - pnorm.eta.2 - pnorm.eta.6 + 2 * pnorm.eta.2 * pnorm.eta.6) *
(- A - B + 2 * (C * D + E * B))
}
pmax(eta, 0) ## gradient cannot be negative.
}
doubleduotrio$linkfun <- function(mu) {
eps <- 1e-10
ok <- mu > 0.25 & mu < 1 - eps
mu[mu <= 0.25] <- 0
mu[mu >= 1 - eps] <- Inf
if(sum(ok)) {
duotriog <- function(d, p) doubleduotrio$linkinv(d) - p
mu[ok] <- sapply(mu[ok], function(mu) {
uniroot(f=duotriog, interval=c(0, 16), p=mu)$root })
}
pmax(mu, 0) ## delta cannot be negative
}
doubleduotrio
}
doublethreeAFC <- function ()
{
doublethreeAFC <- binomial()
doublethreeAFC$link <- "Link for the 3-AFC double test"
doublethreeAFC$linkinv <- function(eta) {
ok <- eta > 0 & eta < 9
eta[eta <= 0] <- 1/9
eta[eta >= 9] <- 1
if(sum(ok)) {
threeAFCg <- function(x, d) dnorm(x - d) * pnorm(x)^2
eta[ok] <- sapply(eta[ok], function(eta) {
(integrate(threeAFCg, -Inf, Inf, d = eta)$value)^2 })
}
pmin(pmax(eta, 1/9), 1) ## restrict to [1/9, 1] - just to be sure
}
doublethreeAFC$mu.eta <- function(eta) {
ok <- eta > 0 & eta < 9
eta[eta <= 0] <- 0
eta[eta >= 9] <- 0
### Note: integration is not reliable for eta > 9
if(sum(ok)) {
threeAFCg <- function(x, d) dnorm(x - d) * pnorm(x)^2
threeAFCgd <- function(x, d)
(x - d) * dnorm(x - d) * pnorm(x)^2
eta[ok] <- sapply(eta[ok], function(eta) {
2 * integrate(threeAFCg, -Inf, Inf, d = eta)$value *
integrate(threeAFCgd, -Inf, Inf, d = eta)$value})
}
pmax(eta, 0) ## gradient cannot be negative.
}
doublethreeAFC$linkfun <- function(mu) {
eps <- 1e-8
ok <- mu > 1/9 & mu < 1 - eps
mu[mu <= 1/9] <- 0
mu[mu >= 1 - eps] <- Inf
if(sum(ok)) {
threeAFCg2 <- function(d, p) doublethreeAFC$linkinv(d) - p
mu[ok] <- sapply(mu[ok], function(mu)
uniroot(threeAFCg2, c(0, 9), p = mu)$root)
}
pmax(mu, 0)
}
doublethreeAFC
}
doubletriangle <- function()
{
doubletriangle <- binomial()
doubletriangle$link <- "Link for the triangle double test"
doubletriangle$linkinv <- function(eta) {
ok <- eta > 0 & eta < 20
eta[eta <= 0] <- 1/9
eta[eta >= 20] <- 1
if(sum(ok))
eta[ok] <-
(pf(q=3, df1=1, df2=1, ncp=eta[ok]^2*2/3, lower.tail=FALSE))^2
pmin(pmax(eta, 1/9), 1) ## restrict to [1/9, 1] - just to be sure
}
# doubletriangle$mu.eta <- function(eta) {
# ok <- eta > 0 & eta < 20
# eta.1 <- eta
# eta.1[eta.1 <= 0] <- 1/9
# eta.1[eta.1 >= 20] <- 1
# if(sum(ok))
# eta.1[ok] <-
# 2*pf(q=3, df1=1, df2=1, ncp=eta.1[ok]^2*2/3, lower.tail=FALSE)
# eta.1 <- pmin(pmax(eta.1, 1/9), 1) ## restrict to [1/9, 1] - just to be sure
# eta[eta <= 0] <- 0
# eta[eta >= 20] <- 0
# if(sum(ok)) {
# d <- eta[ok]
# eta[ok] <- eta.1 * sqrt(2/3) * dnorm(d/sqrt(6)) *
# (pnorm(d/sqrt(2)) - pnorm(-d/sqrt(2)))
# }
#
# pmax(eta, 0) ## gradient cannot be negative.
# }
doubletriangle$mu.eta <- function(eta) {
ok <- eta > 0 & eta < 20
eta[eta <= 0] <- 1/9
eta[eta >= 20] <- 1
if(sum(ok)) {
d <- eta[ok]
eta[ok] <- 2*pf(q=3, df1=1, df2=1, ncp=d^2*2/3, lower.tail=FALSE) *
sqrt(2/3) * dnorm(d/sqrt(6)) *
(pnorm(d/sqrt(2)) - pnorm(-d/sqrt(2)))
}
pmax(eta, 0) ## gradient cannot be negative.
}
doubletriangle$linkfun <- function(mu) {
eps <- 1e-8
ok <- mu > 1/9 & mu < 1 - eps
mu[mu <= 1/9] <- 0
mu[mu >= 1 - eps] <- Inf
if(sum(ok)) {
triangleg2 <- function(d, p) doubletriangle$linkinv(d) - p
mu[ok] <- sapply(mu[ok], function(mu)
uniroot(triangleg2, c(0, 15), p = mu)$root)
}
mu # FIXME: pmax(mu, 0)
}
doubletriangle
}
doubletetrad <- function()
{
doubletetrad <- binomial()
doubletetrad$link <- "Link for the unspecified tetrad test"
doubletetrad$linkinv <- function(eta) {
eps <- 1e-8
ok <- eta > eps & eta < 9
eta[eta <= eps] <- 1/9
eta[eta >= 9] <- 1
if(sum(ok)) {
tetrads.fun <- function(z, delta)
dnorm(z) * (2 * pnorm(z) * pnorm(z - delta) -
pnorm(z - delta)^2)
eta[ok] <- sapply(eta[ok], function(eta) {
(1 - 2*integrate(tetrads.fun, -Inf, Inf, delta=eta)$value)^2 })
}
pmin(pmax(eta, 1/9), 1) ## restrict to [1/9, 1] - just to be sure
}
doubletetrad$mu.eta <- function(eta) {
eps <- 1e-8
ok <- eta > eps & eta < 9
eta[eta <= eps] <- 0
eta[eta >= 9] <- 0
if(sum(ok)) {
Linkinv <- function(eta) {
tetrads.fun <- function(z, delta)
dnorm(z) * (2 * pnorm(z) * pnorm(z - delta) -
pnorm(z - delta)^2)
sapply(eta, function(eta) {
(1 - 2*integrate(tetrads.fun, -Inf, Inf, delta=eta)$value)^2 })
}
eta[ok] <- sapply(eta[ok], function(eta) grad(Linkinv, eta))
### FIXME: Could probably do the integration by hand here.
}
pmax(eta, 0) ## gradient cannot be negative.
}
doubletetrad$linkfun <- function(mu) {
eps <- 1e-8 ## What is the right eps here?
ok <- mu > 1/9 & mu < 1 - eps
mu[mu <= 1/9] <- 0
mu[mu >= 1 - eps] <- Inf
if(sum(ok)) {
doubletetrads <- function(d, p) doubletetrad$linkinv(d) - p
mu[ok] <- sapply(mu[ok], function(mu)
uniroot(doubletetrads, c(0, 9), p = mu)$root)
}
pmax(mu, 0)
}
doubletetrad
}
doubletwoAFC <- function() {
doubletwoAFC <- binomial()
doubletwoAFC$link <- "Link for the 2-AFC double test"
doubletwoAFC$linkinv <- function(eta) {
ok <- eta > 0
eta[!ok] <- 0.25
eta[ok] <- (pnorm(eta[ok] / sqrt(2)))^2
pmin(pmax(eta, 0.25), 1) ## restrict to [0.5, 1] - just to be sure
}
doubletwoAFC$mu.eta <- function(eta) {
eta.linkinv <- eta
ok <- eta > 0
eta[!ok] <- 0
eta.linkinv[!ok] <- 0.25
if(any(ok)) {
sqrt.2 <- sqrt(1/2)
eta.linkinv[ok] <- pnorm(eta.linkinv[ok] / sqrt(2))
eta[ok] <- sqrt(2) * eta.linkinv[ok] * dnorm(eta[ok] * sqrt.2)
}
pmax(eta, 0) ## gradient cannot be negative
}
doubletwoAFC$linkfun <- function(mu) {
ok <- mu > 0.25 & mu < 1
mu[mu <= 0.25] <- 0
mu[mu >= 1] <- Inf
# mu[ok] <- 2 * qnorm(mu[ok])^2
# pmax(mu, 0) ## delta cannot be negative
if(sum(ok)) {
twoafcs <- function(d, p) doubletwoAFC$linkinv(d) - p
mu[ok] <- sapply(mu[ok], function(mu)
uniroot(twoafcs, c(0, 9), p = mu)$root)
}
pmax(mu, 0)
}
doubletwoAFC
}
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.