DirichletProcessMixtureModel/flexloglikelihood.R at master · ashar799/DirichletProcessMixtureModel · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

flexloglikelihood = function(c,Y,mu,S,alpha,That, beta0, betahat, sigma2, lambda2, tau2, K, epsilon, W, beta, ro,D, r, si, Time,N, sig2.dat) {

  disclass <- table(factor(c, levels = 1:K))
  activeclass <- which(disclass!=0)


  Ytemp <- Y
  Ytemp.scaled <- matrix(NA, nrow = N, ncol = D)

  indexrel <- list(0)

  for ( i in 1:length(activeclass)) {
    clust <- which(c == activeclass[i])
    indexrel[[activeclass[i]]] <- which(abs(betahat[activeclass[i],]) > 0.0001)

    if (length(clust)==1){
      Ytemp.scaled[clust,1:D] <- matrix(0, nrow =1, ncol =D)
    } else {
      Ytemp.scaled[clust,1:D] <- scale(Ytemp[clust,1:D], center = TRUE, scale = TRUE)
    }
  }

  loglikelihood <-c(0)

  for (j in 1:length(activeclass)) {

    loglikelihood[j] <- 0
    clust <- which(c == activeclass[j])
    luk1 <- c(0)
    luk2 <- c(0)
    for ( l in 1:length(clust)) {

      relD <- unlist(indexrel[[activeclass[j]]])
      if (length(relD)!= 0){
        Yrelev <- Y[clust[l],relD]
        meanrelev <- mu[activeclass[j],relD]
        precisionrelev <- S[activeclass[j],relD,relD]
      } else {
        Yrelev <- Y[clust[l],1:D]
        meanrelev <- mu[activeclass[j],1:D]
        precisionrelev <- S[activeclass[j],1:D,1:D]
      }

      if (Time[clust[l],2]==1){
        luk1[l] <- log(dMVN(x = as.vector(t(Yrelev)), mean = meanrelev,  Q = precisionrelev))
        luk2[l] <- log(dnorm(x = That[clust[l]], mean = beta0[activeclass[j]] + betahat[activeclass[j],1:D] %*% as.vector(t(Ytemp.scaled[l,1:D])), sd = sqrt(sigma2[activeclass[j]]) ))
      } else{
        luk1[l] <- log(dMVN(x = as.vector(t(Yrelev)), mean = meanrelev,  Q = precisionrelev))
        luk2[l] <- log(dtruncnorm(x = That[clust[l]], a = Time[clust[l],1], b = Inf, mean = beta0[activeclass[j]] + betahat[activeclass[j],1:D] %*% as.vector(t(Ytemp.scaled[l,1:D])), sd = sqrt(sigma2[activeclass[j]]) ))
      }
    }
    loglikelihood[j] <- sum(luk1) + sum(luk2)

  }

  return(sum(loglikelihood))


}