R/M_OEM.R

In DOEM: The Distributed Online Expectation Maximization Algorithms to Solve Parameters of Poisson Mixture Models

#' The M-OEM algorithm replaces M-step with stochastic step, which is used to solve the parameter estimation of Poisson mixture model.
#'
#' @param y is a data vector
#' @param K is the number of Poisson distribution
#' @param alpha0 is the initial value of the mixing weight
#' @param lambda0 is the initial value of the mean
#' @param a represents the power of the reciprocal of the step size
#' @param b indicates that the M-step is not implemented for the first b data points
#'
#' @return M_OEMtime,M_OEMalpha,M_OEMlambda
#' @export
#'

#' @examples
#' library(stats)
#' set.seed(637351)
#' K=5 
#' alpha1=c(rep(1/K,K)) 
#' lambda1=c(1,2,3,4,5) 
#' n=300 
#' U=sample(c(1:n),n,replace=FALSE)
#' y= c(rep(0,n)) 
#' for(i in 1:n){
#' if(U[i]<=0.2*n){
#' y[i] = rpois(1,lambda1[1])} 
#' else if(U[i]>0.2*n & U[i]<=0.4*n){
#' y[i] = rpois(1,lambda1[2])} 
#' else if(U[i]>0.4*n & U[i]<=0.6*n){
#' y[i] = rpois(1,lambda1[3])} 
#' else if(U[i]>0.6*n & U[i]<=0.8*n){
#' y[i] = rpois(1,lambda1[4])}
#' else if(U[i]>0.8*n ){
#' y[i] = rpois(1,lambda1[5])} 
#' }
#' M=5
#' seed=637351
#' set.seed(123) 
#' e=sample(c(1:n),K)
#' alpha0=e/sum(e)
#' lambda0=c(1.5,2.5,3.5,4.5,5.5)
#' a=0.75
#' b=5
#' M_OEM(y,K,alpha0,lambda0,a,b)

M_OEM=function(y,K,alpha0,lambda0,a,b){
n=length(y)
alpha=alpha0
lambda=lambda0
time=system.time(for (t in 0:(n-1)) {
oldalpha=alpha 
oldlambda=lambda
gamma=1/(t+1)^a
w=c(rep(0,K)) 
for (k in 1:K){
w[k]=(oldalpha[k]*exp(-oldlambda[k])* 
oldlambda[k]^ y[t+1])/sum((oldalpha
*exp(-oldlambda)* oldlambda^ y[t+1]))
if(t<b){
alpha[k]=alpha[k]
lambda[k]=lambda[k]
}
else {
alpha[k]= oldalpha [k] +gamma*(w[k]-oldalpha[k])
lambda[k]=oldlambda[k]+ gamma*(w[k]*y[t+1]/oldlambda[k]-w[k])
}
}
cat("alpha",alpha,"\n")
cat("lambda",lambda,"\n")
}
)
return(list(M_OEMtime=time,M_OEMalpha=alpha,M_OEMlambda=lambda))
}