# R/blimMD.R In pks: Probabilistic Knowledge Structures

#### Documented in blimMD

```## Fitting BLIM by minimum discrepancy
blimMD <- function(K, N.R, R = as.binmat(N.R),
betafix = rep(NA, nitems), etafix = rep(NA, nitems),
incrule = c("minimum", "hypblc1", "hypblc2"), m = 1) {

K       <- as.matrix(K)
N.R     <- setNames(as.integer(N.R), names(N.R))  # convert to named int
N       <- sum(N.R)
nitems  <- ncol(K)
npat    <- nrow(R)
nstates <- nrow(K)

## Assigning state K given response R
d.RK <- if(length(which(c(betafix, etafix) == 0))) {
apply(K, 1, function(k) {
RwoK <- t(R) & !k
idx <- which(RwoK, arr.ind=TRUE)
RwoK[idx[idx[, "row"] %in% which(etafix == 0), ]] <- NA

KwoR <- k & !t(R)
idx <- which(KwoR, arr.ind=TRUE)
KwoR[idx[idx[, "row"] %in% which(betafix == 0), ]] <- NA
colSums(RwoK) + colSums(KwoR)
})
} else {
apply(K, 1, function(k) colSums(xor(t(R), k)))
}
d.min <- apply(d.RK, 1, min, na.rm=TRUE)            # minimum discrepancy

i.RK  <- switch(match.arg(incrule),                 # inclusion rule
minimum = (d.RK == d.min) &              !is.na(d.RK),
hypblc1 = replace(1/(1 + d.RK - d.min)^m, is.na(d.RK), 0),
hypblc2 = replace(1/(1 + d.RK)^m,         is.na(d.RK), 0))
m.RK  <- i.RK/rowSums(i.RK) * N.R                   # P(K|R) * N(R)

## Minimum discrepancy distribution
disc.tab <- xtabs(N.R ~ d.min)
disc     <- as.numeric(names(disc.tab)) %*% disc.tab / N

## Distribution of knowledge states
P.K <- colSums(m.RK)/N
names(P.K) <- if(is.null(rownames(K))) as.pattern(K) else rownames(K)

## Careless error and guessing parameters
beta <- eta <- numeric(nitems)
names(beta) <- names(eta) <-
if(is.null(colnames(K))){
make.unique(c("a", letters[(seq_len(nitems) %% 26) + 1])[-(nitems + 1)],
sep="")
}else colnames(K)
for(j in seq_len(nitems)){
beta[j] <- sum(m.RK[which(R[,j] == 0), which(K[,j] == 1)]) /
sum(m.RK[,which(K[,j] == 1)])

eta[j]  <- sum(m.RK[which(R[,j] == 1), which(K[,j] == 0)]) /
sum(m.RK[,which(K[,j] == 0)])
}
beta[is.na(beta)] <- 0
eta[is.na( eta)] <- 0
beta[!is.na(betafix)] <- betafix[!is.na(betafix)]  # reset fixed parameters
eta[!is.na( etafix)] <-  etafix[!is.na( etafix)]

z <- list(discrepancy=c(disc), P.K=P.K, beta=beta, eta=eta,
disc.tab=disc.tab, nstates=nstates, npatterns=npat, ntotal=N,
nerror=NA, method="MD", iter=NA, loglik=NA)
class(z) <- "blim"
z
}
```

## 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.