library(OpenMx)
library(rpf)

## Extended example w/ priors from OpenMx website:

spec <- list()
spec[1:8] <- rpf.drm()  # drm="dichotomous response model"

# replace with your own data
simParam <- sapply(spec, rpf.rparam)
simParam['u',] <- logit(1)     # fix upper bound at 1 for a 3PL equivalent model
data <- rpf.sample(750, spec, simParam)

imat <- mxMatrix(name="item", values=c(1,0,logit(.1),logit(1)),
                free=c(TRUE, TRUE, TRUE, FALSE), nrow=4, ncol=length(spec),
                dimnames=list(names(rpf.rparam(spec[[1]])), colnames(data)))

# label the pseudo-guessing parameters
imat$labels['g',] <- paste('g',1:length(spec),sep="")

# half of the items get a beta prior
betaRange <- 1:(length(spec)/2)

# the other half get a Gaussian prior
gaussRange <- (1+length(spec)/2):length(spec)

# Create matrices that contain only a subset of the parameters from
# the item matrix so the priors are easier to set up.
betaPrior <- mxMatrix(name="betaPrior", nrow=1, ncol=length(betaRange),
               free=TRUE, labels=imat$labels['g',betaRange],
               values=imat$values['g',betaRange])
gaussPrior <- mxMatrix(name="gaussPrior", nrow=1, ncol=length(gaussRange),
               free=TRUE, labels=imat$labels['g',gaussRange],
               values=imat$values['g',gaussRange])


calcBetaParam <- function(mode, strength) {
  a <- mode * strength
  b <- strength - a
  c(a=a, b=b, c=log(beta(a+1,b+1)))
}

guessChance <- c(1/2, 1/3, 1/4, 1/5)
betaParam <- sapply(guessChance, calcBetaParam, strength=5)

# copy our betaParam table into OpenMx row vectors
betaA <- mxMatrix(name="betaA", nrow=1, ncol=length(betaRange), values=betaParam['a',])
betaB <- mxMatrix(name="betaB", nrow=1, ncol=length(betaRange), values=betaParam['b',])
betaC <- mxMatrix(name="betaC", nrow=1, ncol=length(betaRange), values=betaParam['c',])

# implement the math given above
betaFit <- mxAlgebra(2 * sum((betaA + betaB)*log(exp(betaPrior)+1) -
             betaA * betaPrior + betaC), name="betaFit")
betaGrad <- mxAlgebra(2*(betaB-(betaA+betaB) / (exp(betaPrior) + 1)),
             name="betaGrad", dimnames=list(c(),betaPrior$labels))
betaHess <- mxAlgebra(vec2diag(2*exp(betaPrior)*(betaA+betaB) / (exp(betaPrior) + 1)^2), name="betaHess",
             dimnames=list(betaPrior$labels, betaPrior$labels))

# Create a model that will evaluate to the log likelihood of the beta prior
# and provide suitable derivatives for the optimizer.
betaModel <- mxModel(model="betaModel", betaPrior, betaA, betaB, betaC,
             betaFit, betaGrad, betaHess,
             mxFitFunctionAlgebra("betaFit", gradient="betaGrad", hessian="betaHess"))


# These are the prior parameter row vectors.
gaussM <- mxMatrix(name="gaussM", nrow=1, ncol=length(gaussRange), values=logit(guessChance))
gaussSD <- mxMatrix(name="gaussSD", nrow=1, ncol=length(gaussRange), values=.5)

# The single variable Gaussian density and derivatives are well known mathematical results.
gaussFit <- mxAlgebra(sum(log(2*pi) + 2*log(gaussSD) +
             (gaussPrior-gaussM)^2/gaussSD^2), name="gaussFit")
gaussGrad <- mxAlgebra(2*(gaussPrior - gaussM)/gaussSD^2, name="gaussGrad",
             dimnames=list(c(),gaussPrior$labels))
gaussHess <- mxAlgebra(vec2diag(2/gaussSD^2), name="gaussHess",
                   dimnames=list(gaussPrior$labels, gaussPrior$labels))

# Create a model that will evaluate to the log likelihood of the Gaussian prior
# and provide suitable derivatives for the optimizer.
gaussModel <- mxModel(model="gaussModel", gaussPrior, gaussM, gaussSD,
                  gaussFit, gaussGrad, gaussHess,
                  mxFitFunctionAlgebra("gaussFit", gradient="gaussGrad", hessian="gaussHess"))

itemModel <- mxModel(model="itemModel", imat,
           mxData(observed=data, type="raw"),
           mxExpectationBA81(spec),
           mxFitFunctionML())

demo3pl <- mxModel(model="demo3pl", itemModel, betaModel, gaussModel,
                 mxFitFunctionMultigroup(groups=c('itemModel.fitfunction', 'betaModel.fitfunction',
                   'gaussModel.fitfunction')),
                 mxComputeSequence(list(
                   mxComputeEM('itemModel.expectation', 'scores', mxComputeNewtonRaphson()),
                   mxComputeNumericDeriv(),
                   mxComputeHessianQuality(),
                   mxComputeStandardError()
                 )))

demo3pl <- mxRun(demo3pl)
summary(demo3pl)


## Messing around with information matrix and standard errors

demo3pl_xpdse<- mxModel(model="demo3pl_xpdse", itemModel, betaModel, gaussModel,
                 mxFitFunctionMultigroup(groups=c('itemModel.fitfunction', 'betaModel.fitfunction',
                   'gaussModel.fitfunction')),
                 mxComputeSequence(list(
                     mxComputeEM('itemModel.expectation', 'scores', mxComputeNewtonRaphson()),
                       mxComputeOnce('fitfunction', 'information', "sandwich"),
                   ##mxComputeNumericDeriv(),
                   mxComputeHessianQuality(),
                   mxComputeStandardError()
                     )))
demo3pl_xpdse<-mxRun(demo3pl_xpdse)
summary(demo3pl_xpdse)
str(demo3pl_xpdse@output) ## I see standard errors but no information matrix


demo3pl_xpdse2<- mxModel(model="demo3pl_xpdse2", itemModel, betaModel, gaussModel,
                 mxFitFunctionMultigroup(groups=c('itemModel.fitfunction', 'betaModel.fitfunction',
                   'gaussModel.fitfunction')),
                 mxComputeSequence(list(
                     mxComputeEM('itemModel.expectation', 'scores', mxComputeNewtonRaphson()),
                       mxComputeOnce('fitfunction', 'information', "sandwich"),
                   ##mxComputeNumericDeriv(),
                   mxComputeHessianQuality()
                   ##mxComputeStandardError()
                     )))
demo3pl_xpdse2<-mxRun(demo3pl_xpdse2)
summary(demo3pl_xpdse2)
str(demo3pl_xpdse2@output) ## No standard errors or information matrix
## I would assume mxComputeHessianQuality() is based on whatever information matrix is calculated or used?