library(OpenMx)
# specifying the FML function:
simple_FML <- function(model, state) {
  obsCov <- model$data$observed# observed covariance
  nmanif <- ncol(obsCov) # number of observed variables

  Sigma_try_out <- try(mxGetExpected(model,"covariance")) # expected covariance matrix

  if (class(Sigma_try_out)== "try-error"){return(NA)} # if impossible, return NA
  # if it can be built, I am calculating the FML:
  else {
    inv_Sigma_try <- try(solve(Sigma_try_out))
    if (class(inv_Sigma_try)== "try-error"){return(NA)}
    else{
      FML_temp <- log(det(Sigma_try_out))+tr(obsCov%*%inv_Sigma_try)-log(det(obsCov))- nmanif
      if (anyNA(model$data$means)){return(FML_temp)}
      else {
        obsMeans <- model$data$means
        M_try <- try(mxGetExpected(model,"means")) # expected means
        if (class(M_try)== "try-error" ){return(NA)}
        else{FML_temp <- FML_temp + (obsMeans- M_try)%*%inv_Sigma_try%*% t(obsMeans- M_try) #compute FML
        return(FML_temp)}}
  }
  }
}

# Simulating Data:
library(MASS)
set.seed(100)

myData <- mvrnorm (n = 1000, mu= c(0,0,0,0),
                    Sigma= matrix(c(1,0.75,0.75,0,
                                    0.75,1,0.75,0,
                                    0.75,0.75,1,0,
                                    0,0,0,1), byrow = T, 4,4))
colnames(myData) <- c('X1', 'X2', 'X3','X4')


rawCov  <- mxData(observed=cov(myData), type="cov", means = apply(myData, 2,mean), numObs=1000)

##manifest and latent variables:
manifest <- c('X1','X2','X3', 'X4')
latent <- 'LV' 

##Paths
causalPath <- mxPath(from = latent, to = manifest, value = c(0.866,0.866,0.866,0), free = T) 

mancovPath <- mxPath(from = manifest, arrows = 2, free = T,
                     values = c(0.25,0.25,0.25,1))

lvcovPath <- mxPath(from = latent, arrows = 2, value = 1, free = F)


##### calculating base model for comparison #####
baseModel <- mxModel(model = 'Model',
                     type = "RAM", 
                     latentVars = latent, 
                     manifestVars = manifest, 
                     rawCov ,causalPath, mancovPath, lvcovPath,
                     mxPath(from = "one", to = manifest, free = T)
)

fit_baseModel <- mxRun(baseModel)

summary(fit_baseModel)

# the same model but with the ne fitfunction:
myfitfunction <- mxFitFunctionR(simple_FML)
myModel <- mxModel(baseModel, myfitfunction)

fit_myModel <- mxRun(myModel)

summary(fit_myModel)