R/jacobian.R In pks: Probabilistic Knowledge Structures

Documented in jacobian

```jacobian <- function(object, P.K = rep(1/nstates, nstates),
beta = rep(0.1, nitems), eta = rep(0.1, nitems),
betafix = rep(NA, nitems), etafix = rep(NA, nitems))
{
K <- as.matrix(object\$K)
# N <- object\$ntotal
nitems <- ncol(K)
nstates <- nrow(K)
R <- as.matrix(expand.grid(rep(list(0:1), nitems)))
rownames(R) <- as.pattern(R)
npatterns <- nrow(R)

PRKfun <- if(length(which(c(betafix, etafix) == 0))) {
getPRK[["apply"]]
} else {
getPRK[["matmult"]]
}
betanew <- beta
etanew <-  eta
betanew[!is.na(betafix)] <- betafix[!is.na(betafix)]
etanew[!is.na( etafix)] <-  etafix[!is.na( etafix)]
P.R.K <- do.call(PRKfun, list(betanew, etanew, K, R))

# P.R <- as.numeric(P.R.K %*% P.K)
K.star <- K[-nstates, ]
R.star <- R[-npatterns, ]
P.K.star <- P.K[-nstates]
P.R.K.star <- P.R.K[-npatterns, -nstates]
P <- P.R.K.star - outer(P.R.K[-npatterns, nstates], rep(1, nstates - 1))
colnames(P) <- paste("pi", as.pattern(K.star), sep="_")
B <- sapply(1:nitems, function(q) (as.matrix(P.R.K.star[, K.star[, q] == 1]) %*% P.K.star[K.star[, q] == 1] +
P.R.K[-npatterns, nstates] * P.K[nstates]) * (-1/(1 - beta[q]))^R.star[, q] * (1/beta[q])^(1-R.star[, q]))
colnames(B) <- paste("beta", colnames(K), sep="_")
E <- sapply(1:nitems, function(q) (as.matrix(P.R.K.star[, K.star[, q] == 0]) %*% P.K.star[K.star[, q] == 0]) *
(1/eta[q])^R.star[, q] * (-1/(1-eta[q]))^(1-R.star[, q]))
colnames(E) <- paste("eta", colnames(K), sep="_")
J <- cbind(P, B, E)
rownames(J) <- rownames(R.star)
J
}
```

Try the pks package in your browser

Any scripts or data that you put into this service are public.

pks documentation built on July 9, 2023, 5:43 p.m.