##### GG-prior_ln(MM)

gg.pr <- function(data, pattern, dump=F, dump.a=F, dump.z=F, file="mle_gg_pr.txt", file.a="mle_gg_pr_a.txt", file.z="mle_gg_pr_z.txt",
header=null, header2=null) {


# to estimate a_k empirically
a <- apply(data, 1, median)^2 / apply(data, 1, mad)^2
        
# to estimate hyperparameters mu and sigma2 in the lognormal prior using MM
mu <- median(log(a))
sigma2 <- mad(log(a))^2





# pdf of the lognormal prior for a_k
pr <- function(a, mu, sigma2){
  1/a/sqrt(2*pi*sigma2) * exp(-(log(a)-mu)^2/2/sigma2)
}


# log-likelihood for null case (no differential expression)
lp0 <- function(a, a0, v, R, G) {
  n <- ncol(R) + ncol(G)
  value <- lgamma(n*a+a0) - n*lgamma(a) - lgamma(a0) + (a0)*log(v) 
  temp1 <- apply(log(R), 1, sum) + apply(log(G), 1, sum)
  temp2 <- v + apply(R, 1, sum) + apply(G, 1, sum)
  value <- value + (a-1)*temp1 - (n*a+a0)*log(temp2)
  return(value)
}


# log-likelihood for alternative case (differential expression)
lpa <- function(a, a0, v, R, G) {
  nR <- ncol(R)
  nG <- ncol(G)
  n <- nR + nG
  value <- lgamma(nR*a+a0) + lgamma(nG*a+a0) - n*lgamma(a) - 2*lgamma(a0) + (2*a0)*log(v)
  temp1 <- apply(log(R), 1, sum) + apply(log(G), 1, sum)
  temp2 <- v + apply(R, 1, sum)
  temp3 <- v + apply(G, 1, sum)
  value <- value + (a-1)*temp1 - (nR*a+a0)*log(temp2) - (nG*a+a0)*log(temp3)
  return(value)
}


# null loglikelihood
nlogl <- function(aa0v, R, G) {
  value <- sum(lp0(aa0v[1], aa0v[2], aa0v[3], R, G))
  if (is.na(value)) value <- 1e50
  return(-value)
}


# posterior probability of change
z <- function(a, a0, v, p, R, G) {
  p / (p + (1-p)*exp(lp0(a,a0,v,R,G) - lpa(a,a0,v,R,G)))
}


# complete data loglikelihood; maximize wrt a0 and v
clogl_a0v <- function(a0v, R, G, z, p, a) {
  a0 <- a0v[1]
  v <- a0v[2]
  value <- sum(z*lpa(a,a0,v,R,G) + (1-z)*lp0(a,a0,v,R,G))     # + z*log(p) + (1-z)*log(1-p) + (b1-1)*log(p) + (b2-1)*log(1-p)
  if (is.na(value)) value <- 1e50
  return(-value)
}


# gradient of clogl_a0v wrt a0 and v
ga0v <- function(a0v, R, G, z, p, a) {
  a0 <- a0v[1]
  v <- a0v[2]
  nR <- ncol(R)
  nG <- ncol(G)
  n <- nR + nG
  temp1 <- v + apply(R, 1, sum)
  temp2 <- v + apply(G, 1, sum)
  value1 <- sum( z * (digamma(nR*a+a0) + digamma(nG*a+a0) - 2*digamma(a0) + 2*log(v) - log(temp1) - log(temp2)) +
	     (1-z) * (digamma(n*a+a0) - digamma(a0) + log(v) - log(temp1+temp2-v) ) )
  value2 <- sum( z * (2*a0/v - (nR*a+a0) / temp1 - (nG*a+a0) / temp2) + (1-z) * (a0/v - (n*a+a0) / (temp1+temp2-v)) )
  return(c(-value1, -value2))
}
        

# complete data loglikelihood; maximize wrt a_k
clogl_a <- function(a, R, G, z, a0v, p) {
  a0 <- a0v[1]
  v <- a0v[2]
  value <- z*lpa(a,a0,v,R,G) + (1-z)*lp0(a,a0,v,R,G) + log(pr(a, mu, sigma2))
  if (is.na(value)) value <- 1e50
  return(-value)
}


# gradient of clogl_a wrt a_k
ga <- function(a, R, G, z, a0v, p) {
  a0 <- a0v[1]
  v <- a0v[2]
  nR <- ncol(R)
  nG <- ncol(G)
  n <- nR + nG
  value <- -( z * (nR*digamma(nR*a+a0) + nG*digamma(nG*a+a0) - n*digamma(a) + sum(log(R)) + sum(log(G)) - nR*log(v+sum(R)) - nG*log(v+sum(G)))
+ (1-z) * (n*digamma(n*a+a0) - n*digamma(a) + sum(log(R)) + sum(log(G)) - n*log(v+sum(R)+sum(G))) - 1/a *
(1+(log(a)-mu)/sigma2))
  return(value)
} 





n <- nrow(data)
pattern <- strsplit(pattern, " +")[[1]]
con1 <- pattern=="1"
con2 <- pattern=="2"
r <- data[,con1]
g <- data[,con2]


# crude MM estimation of (hyper)parameters
c.sd <- apply(data, 1, sd)
c.mean <- apply(data, 1, mean)
c.cov <- median(c.sd/c.mean, na.rm=T)
c.a <- 1/c.cov^2
c.theta <- c.a/c.mean
c.v <- median(c.theta, na.rm=T) / var( sort(c.theta)[(n%/%4):(n%/%4*3)] )
c.a0 <- c.v*median(c.theta, na.rm=T)
init <- c(c.a, c.a0, c.v)
if (dump) {
  write(header, file, ncolumns=1, append=T)
  write(c("MM", mu, sigma2), file, ncolumns=3, append=T)
  write(c("crude", init, 0.5), file, ncolumns=5, append=T)     # store the crude estimates of (a, a0, v, p)
}
        

# hyperparameters of the beta prior distribution for p
b1 <- 2
b2 <- 2


# to obtain initial estimates of a, a0, v, Z, p using the null loglikelihood
mleinfo <- optim(par=init, fn=nlogl, R=r, G=g, method="L-BFGS-B", lower=c(0.001,0.001,0.001), upper=c(200,200,200))
if (mleinfo$value == -1e50) {
  old <- init
  repeat{
    mleinfo <- optim(par=old, fn=nlogl, R=r, G=g, method="L-BFGS-B", lower=c(0.001,0.001,0.001), upper=c(200,200,200), control=list(maxit=0))
    if (mleinfo$value == -1e50) mleinfo$par <- old
    if (sum(abs( (mleinfo$par-old)/old )<10^-3) == 3) break  else old <- mleinfo$par
  }
}
init2 <- mleinfo$par
if (sum(init2==0.001 | init2==200) > 0) init2 <- init
Z <- z(init2[1], init2[2], init2[3], p=0.5, r, g)
p <- (sum(Z)+b1-1) / (n+b1+b2-2)
init2p <- p
if (dump) write(c("null", init2, init2p), file, ncolumns=5, append=T)     # store the initial estimates of (a, a0, v, p)


# to obtain mles of a
mle_a <- rep(init2[1], n)
for (j in 1:n) {
  mleinfo_a <- optim(par=init2[1], fn=clogl_a, gr=ga, R=t(r[j,]), G=t(g[j,]), z=Z[j], a0v=init2[2:3], p=init2p, method="L-BFGS-B", lower=1e-4, upper=100000)
  if (mleinfo_a$value > -1e50) mle_a[j] <- mleinfo_a$par[1]  else {
      old_a <- init2[1]
      repeat{
        mleinfo_a <- optim(par=old_a, fn=clogl_a, gr=ga, R=t(r[j,]), G=t(g[j,]), z=Z[j], a0v=init2[2:3], p=init2p, method="L-BFGS-B", lower=1e-4, upper=100000, control=list(maxit=0))
        if (mleinfo_a$value == -1e50) mleinfo_a$par[1] <- old_a
        if (abs( (mleinfo_a$par[1]-old_a)/old_a )<10^-3) break  else old_a <- mleinfo_a$par[1]
      }
      mle_a[j] <- mleinfo_a$par[1]
    }
}
if (dump) write(c("a", max(mle_a), min(mle_a), median(mle_a)), file, ncolumns=4, append=T)
if (dump.a) {
  write(header2, file.a, ncolumns=1, append=T)
  write(mle_a, file.a, ncolumns=1, append=T)
  write("", file.a, append=T)
}


# EM algorithm to obtain mles of a0, v and p
mleinfo_a0v <- list()
mleinfo_a0v$par <- init2[2:3]
kk <- 0
repeat{
  old_a0v <- mleinfo_a0v$par
  old_p <- p
  kk <- kk+1
  Z <- z(mle_a, old_a0v[1], old_a0v[2], old_p, r, g)
  p <- (sum(Z)+b1-1) / (n+b1+b2-2)     # MLE of p
  mleinfo_a0v <- optim(par=old_a0v, fn=clogl_a0v, gr=ga0v, R=r, G=g, z=Z, p=old_p, a=mle_a, method="L-BFGS-B", lower=c(0.001,0.001), upper=c(200,200))
  if (mleinfo_a0v$value == -1e50) mleinfo_a0v$par <- old_a0v
  if (dump) write(c(kk, mleinfo_a0v$value, mleinfo_a0v$par, p), file, ncolumns=5, append=T)
  if (sum( abs((mleinfo_a0v$par-old_a0v)/old_a0v)<10^(-3) ) + sum( abs((p-old_p)/old_p)<10^(-3) ) ==3) break
  if (kk==100) break
}
if (dump) write("", file, append=T)
mle_a0v <- mleinfo_a0v$par
mle_p <- p
Z <- z(mle_a, mle_a0v[1], mle_a0v[2], mle_p, r, g)
if (dump.z) {
  write(header2, file.z, ncolumns=1, append=T)
  write(Z, file.z, ncolumns=1, append=T)
  write("", file.z, append=T)
}

list(mu=mu, sigma2=sigma2, a=mle_a, a0=mle_a0v[1], v=mle_a0v[2], p=mle_p, z=Z)

}