Nothing
lnl.ordinal <- function(param, y, X, id, model, link, rn, start.sigma = FALSE){
if (start.sigma) gradient <- hessian <- TRUE
# F, f, e la fonction de distribution de proba et ses derivees
# premiere et seconde
if (link == "probit"){
F <- pnorm
f <- dnorm
e <- function(x) - x * f(x)
}
if (link == "logit"){
F <- function(x) exp(x)/(1+exp(x))
f <- function(x) exp(x)/(1+exp(x))^2
e <- function(x) (exp(x) - exp(3*x))/(1+exp(x))^4
}
K <- ncol(X)
n <- length(unique(id))
J <- length(unique(y))
beta <- param[1L:K]
mu <- c(-100, 0, param[(K+1L):(K+J-2L)], 100)
bX <- as.numeric(crossprod(t(X), beta))
if (model == "pooling"){
z1 <- mu[y+1] - bX ; z2 <- mu[y] - bX
F1 <- F(z1) ; F2 <- F(z2)
Li <- F1 - F2
}
if (model == "random"){
sigma <- param[K+J-1L]
Pitr <- lapply(rn$nodes, function(x)
F(mu[y+1] - bX - sigma * x) - F(mu[y] - bX - sigma * x))
Pir <- lapply(Pitr, function(x) tapply(x, id, prod))
Li <- Reduce("+", mapply("*", Pir, rn$weights, SIMPLIFY = FALSE))/sqrt(pi)
}
lnL <- sum(log(Li))
W <- matrix(0, nrow(X), J)
# The first two cols are not estimated coefficients (-infty and 0)
# so remove them
W1 <- (col(W) == (y + 1))[, -c(1, 2)] ; W2 <- (col(W) == y)[, -c(1, 2)]
if (model == "pooling"){
f1 <- f(mu[y+1] - bX) ; f2 <- f(mu[y] - bX)
gmu <- (W1 * f1 - W2 * f2)
gb <- - (f1 - f2) * X
gradi <- cbind(gb, gmu) / Li
}
if (model == "random"){
gitr <- mapply(
function(w, x, p){
f1 <- f(mu[y+1] - bX - sigma * x) ; f2 <- f(mu[y] - bX - sigma * x)
F1 <- F(mu[y+1] - bX - sigma * x) ; F2 <- F(mu[y] - bX - sigma * x)
w * as.numeric(p[as.character(id)]) *
(f1 * cbind(- 1, W1, - x) -
f2 * cbind(- 1, W2, - x)) / (F1 - F2)
},
rn$weights, rn$nodes, Pir, SIMPLIFY = FALSE)
g <- Reduce("+", gitr)
gradi <- cbind(X * g[, 1], g[, -c(1)])/as.numeric(Li[as.character(id)])/sqrt(pi)
}
if (model == "pooling"){
e1 <- e(mu[y+1] - bX) ; e2 <- e(mu[y] - bX)
M1 <- cbind(- X, W1)
M2 <- cbind(- X, W2)
H <- crossprod(M1 * e1 / Li, M1) - crossprod(M2 * e2 / Li, M2) - crossprod(gradi)
}
if (model == "random"){
H <- mapply(
function(w, x, p){
p <- p / Li
P <- p[as.character(id)]
p <- as.numeric(p)
P <- as.numeric(P)
z1 <- mu[y+1] - bX - sigma * x
z2 <- mu[y] - bX - sigma * x
f1 <- f(z1) ; f2 <- f(z2)
e1 <- e(z1) ; e2 <- e(z2)
F1 <- F(z1) ; F2 <- F(z2)
M1 <- cbind(- X, W1, - x) ; M2 <- cbind(- X, W2, - x)
git <- (f1 * M1 - f2 * M2) / (F1 - F2)
gi <- apply(git, 2, tapply, id, sum)
H1 <- crossprod(p * gi, gi)
H2 <- crossprod(P * (e1 * M1) / (F1 - F2), M1) -
crossprod(P * (e2 * M2) / (F1 - F2), M2) -
crossprod( git * sqrt(P))
(H1 + H2) * w},
rn$weights, rn$nodes, Pir, SIMPLIFY = FALSE)
H <- Reduce("+", H) / sqrt(pi) - crossprod(apply(gradi, 2, tapply, id, sum))
}
if (start.sigma){
gradit <- (- f1 + f2 ) / Li
hessit <- (e1 - e2) / Li - gradit^2
mu <- tapply(gradit, id, sum) / tapply(hessit, id, sum)
return(sqrt(2) * sd(mu))
}
attr(lnL, "gradient") <- gradi
attr(lnL, "hessian") <- H
lnL
}
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.