spfit <- function(x,y, ngrid=200, nrecurse=50)
 {
 ########################################################################
 ##
 ##  spfit: R code to fit a semiparametric hierarchical mixture model 
 ##         for differential expression analysis.
 ##         Tested on R Version 1.7.1
 ##
 ##  Version 1.1  November, 2003.
 ##
 ##  Copyright (C) 2003 Michael A. Newton  (newton@stat.wisc.edu)
 ##
 ## Input:
 ##
 ##  x       matrix of **RAW** expression levels, group 1 
 ##          (rows = genes and columns equal replicate arrays)
 ##  y       matrix of same for group 2
 ##
 ##          It is assumed that x and y represent data after any
 ##          background adjustment and normalization. They
 ##          are on the raw rather than log scale.
 ##
 ##  ngrid   size of even grid on scale of theta = log(expected expression)
 ##  nrecurse  number of EM iterations
 ##
 ##
 ## Returns a list:
 ##
 ##  post     ngene x 3 matrix holding posterior probs of each hypothesis
 ##           for each gene; hypotheses: mu1=mu2, mu1<mu2, mu1>mu2.
 ##
 ##  grid     a vector holding grid points for evaluation of base 
 ##           distribution pi which governs log(mu)'s as random effects
 ##
 ##  pi       probability vector on grid
 ##
 ##  pmix     estimated overall proportions of mu1=mu2, mu1<mu2, mu1>mu2
 ##
 ##  shapes   estimated shape parameters of Gamma observation component
 ##
 ##
 ##  While running, the code prints the EM recursion count and the
 ##  current value of pmix, which should stabilize over iterations.
 ##  (To do: The internal quantity Del may be manipulated to derive another
 ##  convergence diagnostic.)
 ##
 ##  WARNING: As with any EM algorithm, this one may be slow 
 ##     (order ngene x ngrid x nrecurse).
 ##
 ##
 ##  Model:
 ##
 ##   x_i, y_i conditionally independent Gamma distributed measurements
 ##   with condition-specific shape and with gene/condition specific mean.
 ##   The means (mu1,mu2) of x_i and y_i are random effects from the
 ##   following mixture.  A three sided die is tossed with probabilities
 ##   in pmix.  If 0, then mu1=mu2 and the common value (logged) comes
 ##   from pi.  If 1, then mu1<mu2 and the values come as order statistics
 ##   from two independent draws from pi.  If 2, mu1>mu2 and the values
 ##   similarly come as order stats. from two independent draws from pi.
 ##
 ##  Method: EM on grid of log(mu) values to estimate the base dist. 
 ##  pi and the mixing proportions pmix.  Also returns gene-specific 
 ##  posterior distribution over three patterns of expression.
 ##
 ##  This version is motivated by Technical Report 1074, Jan 03, UW Madison
 ##  Department of Statistics, but it is a major  revision (July 03) of that
 ##  method in that it works with a 1-d mixing distribution. 
 ##  
 ##  Send any/all inquiries to newton@stat.wisc.edu
 ##
 ############################################################################
 ##
 ## This program is free software; you can redistribute it and/or modify it
 ## under the terms of the GNU General Public License as published by the
 ## Free Software Foundation; either version 2, or (at your option) any
 ## later version.
 ##
 ## These functions are distributed in the hope that they will be useful,
 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 ## GNU General Public License for more details.
 ##
 ## The text of the GNU General Public License, version 2, is available
 ## as http://www.gnu.org/copyleft or by writing to the Free Software
 ## Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 ##
 #########################################################################


 ngene <- nrow(x)
 nrep1 <- ncol(x)
 nrep2 <- ncol(y)

 cv1 <- ( apply(x,1,cv) )^2     
 cv2 <- ( apply(y,1,cv) )^2    
 
 a1 <- 1/( mean(cv1) )  # moment estimate of shape parameter for x's
 a2 <- 1/( mean(cv2) )  # moment estimate of shape parameter for y's

 # initiate the mixing distribution on a grid
 tmp <- range( c( log(x), log(y) ) )   
 grid <- seq( tmp[1], tmp[2], length=ngrid )  
 delta <- grid[2] - grid[1]
 pi <- rep(1/ngrid, ngrid)  # initial estimate of base mixing distribution
 pmix <- c(1/3,1/3,1/3)

 one1 <- matrix(1,nrep1,1)
 one2 <- matrix(1,nrep2,1)

 # Joint prob( data_j | mixing parameter ), across whole grid
 sumx <- drop( x %*% one1 )
 sumlx <- drop( log(x) %*% one1 )
 llikx <-  t( nrep1*( a1*( log(a1) - grid ) - lgamma(a1) ) + 
              t( (a1-1) * sumlx - sumx %o% ( a1 * exp(-grid ) ) ) )
 sumy <- drop( y %*% one2 )
 sumly <- drop( log(y) %*% one2 )
 lliky <-  t( nrep2*( a2*( log(a2) - grid ) - lgamma(a2) ) + 
              t( (a2-1) * sumly - sumy %o% ( a2 * exp(-grid ) ) ) )

 #  Matrix to hold posterior probabilities of each pattern

 post <- matrix(NA,ngene,3) 
 
 # Recursion to get NPMLE of mixing distribution

 pi.old <- pi
 pred <- numeric(3)

 for( i in 1:nrecurse )
  {
   # loop through genes
   pi.hat <- rep(0,ngrid) 
   for( j in 1:ngene )
    {
     # 1-d marginal posterior quantities

     postx <- llikx[j,] + log(pi)
     cx <- max(postx)
     postx <- exp(postx-cx)
     sx <- sum( postx )
     margx <- cx + log( sx )
     postx <- postx/sx
     cdfx <- cumsum(postx)

     posty <- lliky[j,] + log(pi)
     cy <- max(posty)
     posty <- exp(posty-cy)
     sy <- sum( posty )
     margy <- cy + log( sy )
     posty <- posty/sy
     cdfy <- cumsum(posty)

     postxy <- llikx[j,] + lliky[j,] + log(pi)
     cxy <- max(postxy)
     postxy <- exp(postxy-cxy)
     sxy <- sum( postxy )
     margxy <- cxy + log( sxy )
     postxy <- postxy/sxy

     # (log) marginal predictive prob. masses under different hypotheses

     pred[1] <- margxy
     pred[2] <- log(2) + margx + margy + log( max(0,sum( postx*(1-cdfy)  )) )
     pred[3] <- log(2) + margx + margy + log( max(0,sum( posty*(1-cdfx)  )) )
     
     tmp <- max(pred)
     marg <- tmp + log( sum( pmix*exp(pred-tmp) ) )   # fully marginal prob.
 

     part0 <-  pmix[1]*exp( margxy - marg )*( postxy + pi )

     part1 <-  2*pmix[2]*exp(margx+margy-marg)*( posty*cdfx + postx*(1-cdfy) )
 
     part2 <-  2*pmix[3]*exp(margx+margy-marg)*( postx*cdfy + posty*(1-cdfx) )

     # psi = Prob[ U_i = mu | x_i,y_i ] + Prob[ V_i = mu | x_i,y_i ] 

     psi <- part0 + part1 + part2

     # Update estimate of pi using data on gene j

     pi.hat <- pi.hat + psi

     # Record posterior probability info

     tmp <- max(pred)
     postz <-  pmix*exp(pred-tmp) 
     postz <- postz/sum(postz)
     post[j,] <- postz
 
     }

    # rescale estimated mixing dist.

    pi.hat <- pi.hat/sum(pi.hat)
    pi.old <- pi
    pi <- pi.hat
    
    pmix <- apply(post,2,mean)   # E step on mixing proportions
 
    print( c( i, round( pmix, 3 ) ) )
   
    Del <- ngene*( pi/pi.old - 1 ) # Could be used for conv. diagnostic
    pi.old <- pi
  }

 # label stuff
 gene.labs  <- dimnames(x)[[1]]
 hyp.labs1 <-  c("P(mu1=mu2|x,y)", "P(mu1<mu2|x,y)", "P(mu1>mu2|x,y)" )
 hyp.labs2 <-  c("P(mu1=mu2)", "P(mu1<mu2)", "P(mu1>mu2)" )

 dimnames( post ) <- list( gene.labs, hyp.labs1 )
 names(pmix) <- hyp.labs2

 res <- list( post=post, pi=pi, pmix=pmix, grid=grid, shapes=c(a1,a2) )
 return(res)
 }

 cv <- function(x) { return(sqrt(var(x))/mean(x)) }