LIS_ML <- function(data,K){
  
  #log likelihood function
  ll_lis <- function(arg){
    pi_raw <- arg[1:(K-1)]
    alpha <- arg[K:(K+P-1)]
    beta <- arg[(K+P):(K+2*P-1)]
    
    temp = sum(exp(pi_raw))+exp(-5)
    pi = c(exp(pi_raw),exp(-5))/temp
    
    mu <- t(seq(1:K))%*%pi
    sigma <- (t(seq(1:K)-mu)^2%*%pi)^0.5
    theta <- (seq(1:K)-mu)/sigma
    
    log_lik <- function(x,theta,alpha,beta){
      ps <- function(x){
        bernoulli_lik <- function(theta){
          prob <- 1/(1+exp(-1.7*alpha*(theta-beta)))
          temp <- function(x,prob) log(dbinom(x,1,prob))
          return(sum(mapply(temp,x,prob)))
        }
        return(sapply(theta,bernoulli_lik))
      }
      return(sapply(x,ps))
    }
    
    temp <- log_lik(x,theta,alpha,beta)
    ll <- sum(log(apply(exp(log(pi)+temp),2,sum)))
    return(ll)
  }
  
  #data
  N <- nrow(data)
  P <- ncol(data)
  x <- as.list(data.frame(t(as.matrix(data))))
  
  #initial values
  arg <- c(rep(0,K-1),rep(1,P),rep(0,P))
  
  #optimization
  opt <- optim(arg,ll_lis,control=list(fnscale=-1),method="BFGS",hessian = TRUE)
  
  #parameter output
  pi_raw <- opt$par[1:(K-1)]
  temp <- c(pi_raw,-5)
  pi <- exp(temp)/(sum(exp(temp)))
  alpha <- opt$par[K:(K+P-1)]
  beta <- opt$par[(K+P):(K+2*P-1)]
  
  #standard error output
  se <- (-diag(solve(opt$hessian)))^0.5
  pi.se <- se[1:(K-1)]
  alpha.se <- se[K:(K+P-1)]
  beta.se <- se[(K+P):(K+2*P-1)]
  
  return(list(pi=pi,alpha=alpha,beta=beta,pi.se=pi.se,alpha.se=alpha.se,beta.se=beta.se,opt=opt))
}

LIS_score <- function(data,result){
  x <- as.list(data.frame(t(as.matrix(data))))
  alpha <- result$alpha
  beta <- result$beta
  pi2 <- result$pi
  K <- length(pi2)
  mu <- t(seq(1:K))%*%pi2
  sigma <- (t(seq(1:K)-mu)^2%*%pi2)^0.5
  theta <- (seq(1:K)-mu)/sigma
  
  log_lik <- function(x,theta,alpha,beta){
    ps <- function(x){
      bernoulli_lik <- function(theta){
        prob <- 1/(1+exp(-1.7*alpha*(theta-beta)))
        temp <- function(x,prob) log(dbinom(x,1,prob))
        return(sum(mapply(temp,x,prob)))
      }
      return(sapply(theta,bernoulli_lik))
    }
    return(sapply(x,ps))
  }
  
  loglik <- log_lik(x,theta,alpha,beta)
  theta_ind <- t(exp(log(pi2)+loglik))/(apply(exp(log(pi2)+loglik),2,sum))
  score <- apply(theta_ind,1,which.max)
  score <- as.vector(score)
  return(score)
}