R/FGSN.EM.R

In a-mahdavi/SSMFGSN.EM: EM algorithms for the SSMFGSN distributions

FGSN.EM <- function(y, xi, s, la, iter.max=200, tol=10^-6){
  m <- length(la) ; a=seq(1,2*m-1,by=2)
  n <- length(y)
  dif <- 1
  count <- 0
  alp <- as.matrix(la/s^a)
  y.xi <- y-xi
  a.y.xi <- t(outer(y.xi,a,'^'))
  alp.a.y.xi <- as.numeric(t(alp)%*%a.y.xi)
  while ((dif > tol) && (count <= iter.max)) {
    # E step
    s3 <- alp.a.y.xi + dnorm(alp.a.y.xi)/pnorm(alp.a.y.xi)
    LL <- sum(log(2/s*dnorm(y.xi/s)*pnorm(alp.a.y.xi))) # log-likelihood function
    # MCE steps
    xi <-  optim(par=xi,function(x){
      y.xi <- y-x
      a.y.xi <- t(outer(y.xi,a,'^'))
      alp.a.y.xi <- as.numeric(t(alp)%*%a.y.xi)
      return(-sum(log(2/s*dnorm(y.xi/s)*pnorm(alp.a.y.xi))))
    },method= "BFGS")$par
    y.xi <- y-xi
    s <- sqrt(1/n*sum(y.xi^2))
    a.y.xi <- t(outer(y.xi,a,'^'))
    s3.a.y.xi <- matrix(0,m,n)
    for(i in 1:n)
      s3.a.y.xi[,i] <- s3[i]*a.y.xi[,i]
    alp <- solve(a.y.xi%*%t(a.y.xi))%*%rowSums(s3.a.y.xi)
    alp.a.y.xi <- as.numeric(t(alp)%*%a.y.xi)
    LL.new <- sum(log(2/s*dnorm(y.xi/s)*pnorm(alp.a.y.xi))) # log-likelihood function
    count <- count +1
    dif <- abs(LL.new/LL-1)
  }
  aic <- -2 * LL.new + 2 * (2+m)
  bic <- -2 * LL.new + log(n) * (2+m)
  list(xi=xi, sigma=s, la=alp*s^a, loglik=LL.new, AIC=aic, BIC=bic, iter=count)
}