R/dmix.R

In fitmix: Finite Mixture Model Fitting of Lifespan Datasets

#' The mixture distribution
#' @description Computing probability density function for the well-known mixture models.
#' @param lifespan Vector of samples
#' @param model choice of one of the mixture models; \code{gompertz}, \code{log-logistics}, \code{log-normal}, and \code{weibull}.
#' @param K number of components
#' @param param Vector of weight \eqn{\omega}, shape \eqn{\alpha}, and scale \eqn{\beta} parameters.
#' @return A vector of the same length as lifespan data, given the pdf of the one of the mixture models computed at \code{lifespan}.
#' @export
#' @importFrom stats dlnorm dweibull kmeans median nlminb optim plnorm pweibull qlnorm rmultinom runif rweibull
#' @examples lifespan<-seq(0,30,0.2)
#' K<-2
#' weight<-c(0.6,0.4)
#' alpha<-c(0.5,1)
#' beta<-c(1,0.5)
#' param<-c(weight,alpha,beta)
#' dmix(lifespan, "log-logistic", K, param)

dmix<-function(lifespan,model,K,param){
  if(model!="gompertz" & model!="log-logistic" & model!="log-normal" & model!="weibull"){
    stop ("No model detected!")}
  y<-rep(NA,length(lifespan));

  if(model=="gompertz"){
    func=function(x,par){a=par[1]; b=par[2]; b*exp(a*x)*exp(-(exp(a*x)-1)*b/a)}
  }
  if(model=="log-logistic"){
    func=function(x,par){a=par[1]; b=par[2]; (a*b^(-a)*x^(a-1))/(((x/b)^a +1)^2)}
  }
  if(model=="log-normal"){
    func=function(x,par){a=par[1]; b=par[2]; dlnorm(x,a,sdlog=b)}
  }

  if(model=="weibull"){
    func=function(x,par){a=par[1]; b=par[2]; dweibull(x,a,scale=b)}
  }
  kappa<-param[1:K]
  alpha<-param[(K+1):(2*K)]
  beta<-param[(2*K+1):(3*K)]
  lambda<-param[(3*K+1):(4*K)]
  for (i in 1:length(lifespan)){
    pdf<-0
    for (j in 1:K){
      pdf<-pdf+kappa[j]*func(lifespan[i],c(alpha[j],beta[j],lambda[j]))
    }
    y[i]<-pdf
  }
  return(y)
}