DICdat=function(n=1000,nburn=100,k=2,data=scan("Galaxy.data")){
 #
 # Plain Gibbs estimation of the Gaussian mixture model
 #
 # Tidying things up in Glasgow on 02/11/03_______________
 #	adding overall density estimate
 #	adding estimated overall MAP estimate of (z,theta)
 # 	adding estimated marginal MAP estimate of theta
 # 	adding new DIC's and removing old DIC's
 #
 # Completing in Paris and Grenoble_________________
 # addition of DIC5 09/25/02 [removed on 02/12/03]
 # addition of DIC2 09/26/02 [removed on 02/13/03]
 # addition of a burn in stage with nburn steps
 # estimated log-variances
 #
 # addition of DICcomp and DEVcomp 04/10/02
 # correction d'une erreur dans order order le 11/02/03
 #
 # Started on Tenerife, 06/07/02______________

 # Prior distributions 
 #  Dirichlet(alpha) on the weights (p1,...,pk)
   alpha=rep(1,k)
 #  Nnormal N(xi,sigma^2/enne) on the mu's
   xi=rep(0,k)
   enne=rep(0.1,k)
 #  Inverse gamma Ga(nu/2,esse2/2) on the sigma^2's
 # the nu's must be greater than 2
   nu=rep(1.,k)
   esse2=rep(1,k)
 
 # Data
  obs=data
  nobs=length(obs)

 # EM algorithm for initialisation purposes only!
  emresu=EMalgo(obs,k,25)
  
 #  DIC criteria
  DIC1=0
  DICobs=0
  DICAz=0
  DICAt=0
  DICB=0
  DICE=0

 # DICA(t and z) corresponds to the posterior mean deviance
 # when the complete likelihood is used to define the deviance
 # DICB corresponds to the deviance taken at the posterior mean theta
 # Both are related to DIC=c in Flo's manuscript of Feb, 2003.....

 # Initialisation
  uno=t(rep(1,nobs))	# vector of 1's of size nobs
  dos=rep(1,k) 		# vector of 1's of size k
  xranJ=seq(min(obs)-0.5,max(obs)+0.5,(1+diff(range(obs)))/1000)
  nranJ=length(xranJ)
  trio=t(rep(1,nranJ))	# vector of 1's of size nranJ

 # normal normalising constant
  demlog2pi=0.5*nobs*log(2*pi)

 # More initialisation for empirical joint MAP
  zbest=0*(1:nobs)
  pebest=rep(1,k)/k
  mubest=rep(0,k)
  sibest=rep(1,k)
  best=-10000000

 # And more for empirical marginal MAP
  pebesM=emresu$p
  mubesM=emresu$mu
  sibesM=emresu$sig
  compas=2*emresu$like+lprior(emresu$p,emresu$mu,emresu$sig,alpha,esse2,enne,nu,xi)
  besM=-1000000

 # First values of the parameters
  p=rexp(k)
  p=p/sum(p)
  mu=mean(obs)+rnorm(k,0,0.5*diff(range(obs)))
  sig=var(obs)/rgamma(k,2.,2.)

 # First allocation vector
  prob=(p%*%uno)*exp(-0.5*(dos%*%t(obs)-mu%*%uno)^2/sig)/sqrt(sig%*%uno)
  prob=t(t(prob)/apply(prob,2,sum))
  z=(k+1)-apply((runif(nobs)<apply(prob,2,cumsum)),2,sum)

 # More initialisation
  avep=0*p
  avmu=0*mu
  aves=avmu
  avez=0*obs
  avmz=0*t(z%*%t(dos))
  avdens=0*(1:nobs)
  plodens=0*(1:nranJ)

 # Burn in stage
  for (i in 1:nburn){

     fx=apply(dos%*%t(z)==(1:k),1,sum)			     # cpt sizes
     sx=apply( ( t( t(z%*%t(dos))==(1:k) ) )*obs,2,sum)      # cpt sums
     ss=apply( ( t( t(z%*%t(dos))==(1:k) ) )*(obs)^2 ,2,sum) # cpt s of s

     p=rgamma(k,alpha+fx,1)
     p=p/sum(p)
     mu=sqrt(sig/(enne+fx))*rnorm(k)+sx/(enne+fx)
     sig=0.5*(esse2+ss-sx^2/(enne+fx)+enne*xi*(fx*xi-2*sx)/(enne+fx))/rgamma(k,nu/2+0.5*fx)

  # Ordering according to mu's
    p=p[order(mu)]
    sig=sig[order(mu)]
    mu=sort(mu)

  # New allocations
    prob=((dos%*%t(obs))-(mu%*%uno))^2/(sig%*%uno)
    prob=(p%*%uno)*exp(-0.5*prob)/sqrt(sig%*%uno)
    prob=t(t(prob)/apply(prob,2,sum))
    z=(k+1)-apply((runif(nobs)<apply(prob,2,cumsum)),2,sum)

  # Complete likelihood+prior computation
    comp9=2*log(p[z])-(obs-mu[z])^2/sig[z]-log(sig[z])
    like=sum(comp9)+lprior(p,mu,sig,alpha,esse2,enne,nu,xi)

  # (overall) MAP actualisation
    if (best<like){
	best=like
	zbest=z
	mubest=mu
	pebest=p
	sibest=sig
	}

  # Marginal likelihood+prior computation
    prob=apply((p%*%uno)*exp(-0.5*(dos%*%t(obs)-mu%*%uno)^2/sig)/sqrt(sig%*%uno),2,sum)
    like=2*sum(log(prob))+lprior(p,mu,sig,alpha,esse2,enne,nu,xi)

  # (marginal) MAP actualisation
    if (besM<like){
	besM=like
	mubesM=mu
	pebesM=p
	sibesM=sig
	}
    }

 # Main Gibbs loop
  for (i in 1:n){

  # -----------------Gibbs part----------------------------------------

     fx=apply(dos%*%t(z)==(1:k),1,sum)			     # cpt sizes
     sx=apply( ( t( t(z%*%t(dos))==(1:k) ) )*obs,2,sum)      # cpt sums
     ss=apply( ( t( t(z%*%t(dos))==(1:k) ) )*(obs)^2 ,2,sum) # cpt s of s

   # Simulating the parameters
     p=rgamma(k,alpha+fx,1)
     p=p/sum(p)
     mu=sqrt(sig/(enne+fx))*rnorm(k)+sx/(enne+fx)
     sig=0.5*(esse2+ss-sx^2/(enne+fx)+enne*xi*(fx*xi-2*sx)/
     	(enne+fx))/rgamma(k,nu/2+0.5*fx)

  # Ordering according to mu's (usually the best constraint)
    p=p[order(mu)]
    sig=sig[order(mu)]
    mu=sort(mu)

  # New allocation
    prob=((dos%*%t(obs))-(mu%*%uno))^2/(sig%*%uno)
    prob=(p%*%uno)*exp(-0.5*prob)/sqrt(sig%*%uno)
    prob=t(t(prob)/apply(prob,2,sum))

    z=(k+1)-apply((runif(nobs)<apply(prob,2,cumsum)),2,sum)
    matz=(t(z%*%t(dos))==(1:k))
  
  # -----------------Averages and MAPs---------------------------------------

    avep=avep+p
    avmu=avmu+mu
    aves=aves+log(sig)
    avez=avez+z
    avmz=avmz+matz
  
  # Complete likelihood computation
    comp9=2*log(p[z])-(obs-mu[z])^2/sig[z]-log(sig[z])
    like=sum(comp9)+lprior(p,mu,sig,alpha,esse2,enne,nu,xi)

  # (overall) MAP actualisation
    if (best<like){
	best=like
	zbest=z
	mubest=mu
	pebest=p
	sibest=sig
	}

  # Marginal likelihood+prior computation
    prob=apply((p%*%uno)*exp(-0.5*(dos%*%t(obs)-mu%*%uno)^2/sig)/sqrt(sig%*%uno),2,sum)
    like=2*sum(log(prob))+lprior(p,mu,sig,alpha,esse2,enne,nu,xi)

  # (marginal) MAP actualisation
    if (besM<like){
	besM=like
	mubesM=mu
	pebesM=p
	sibesM=sig
	}

  # -----------------DIC part------------------------------------------------

  # To compute DIC=c
  # Estimating "\bar{theta}(y, z^{(l)})" defined in the manuscript page 5
  # The estimations will be in pt, mut, sigt, 
  # to compute the second term in DICcomp, ie DICB

  # mt corresponds to m=j^{(l)} page 5, the number of observations
  # allocated to class j
   mt=apply(matz,1,sum)

  # The following operation is to avoid dividing by 0 when computing muhat
  # and s2t that correspond to \hat{\mu}=j and \hat{s}^2=j.
  # Note that if mt=0 then matz...=0 so taking mt=1 does not change the result
   mt[mt==0]=1

  # Computing \hat{\mu}=j^{(l)} and \hat{s}=j^2{(l)}
   muhat=(matz%*%obs)/mt
   s2t=(matz%*%((obs-muhat[z])^2))

  # formulas page 5 using the parameters of the prior
  # WARNING : I removed the earlier "-2" 
   sigt=(esse2+s2t+ enne*mt*(muhat-xi)^2/(enne+mt))/(nu+mt)
   mut=(enne*xi+mt*muhat)/(enne+mt)
   pt=(mt+alpha)/(sum(alpha)+nobs)

   tinorm=0*(1:nobs)
   compt=0*(1:nobs)
   for (j in 1:k){
       temp=p[j]*exp(-(obs-mu[j])^2/(2*sig[j]))/sqrt(sig[j])
       logtemp=log(p[j])-((obs-mu[j])^2/(2*sig[j]))-0.5*log(sig[j])
       tinorm=tinorm+temp
       compt=compt+ temp * logtemp
       }
   compt=compt/tinorm
   DICAt=DICAt+sum(compt)-demlog2pi

  # Osolete________________________________________________=
  # Computing DICAz (equation (11) page 4) and DICB
  # logcomp=log(p[z])-(obs-mu[z])^2/(2*sig[z])-.5*log(sig[z])
  # DICAz=DICAz+sum(logcomp)-demlog2pi

   logcomp=-(obs-mut[z])^2/(2*sigt[z])-.5*log(sigt[z])
   DICE=DICE+sum(logcomp)-demlog2pi
   logcomp=logcomp+log(pt[z])
   DICB=DICB+sum(logcomp)-demlog2pi

  # Likelihood conditional on the z's, [DIC | z]
   logcomp=-((obs-mu[z])^2/sig[z])-log(sig[z])
   DIC1=DIC1+0.5*sum(logcomp)-demlog2pi

  # Computing the DIC2 criterion of [1], observed likelihood
   prob=apply((p%*%uno)*exp(-0.5*(dos%*%t(obs)-mu%*%uno)^2/sig)/sqrt(sig%*%uno),2,sum) 
   DICobs=DICobs+sum(log(prob))-demlog2pi

   avdens=avdens+prob/sqrt(2*pi)

   prob=apply((p%*%trio)*exp(-0.5*(dos%*%t(xranJ)-mu%*%trio)^2/sig)/sqrt(sig%*%trio),2,sum)
   plodens=plodens+prob/sqrt(2*pi)

   } # end of the main loop

   # averages
   avep=avep/n
   avmu=avmu/n
   aves=aves/n
   aves=exp(aves)
   avez=avez/n
   avmz=avmz/n
   avdens=avdens/n
   plodens=plodens/n

   # scaling (pre-)DICs
   DIC1=-2*DIC1/n
   DICobs=-2*DICobs/n
   #DICAz=-2*DICAz/n
   DICAt=-2*DICAt/n
   DICB=-2*DICB/n
   DICE=-2*DICE/n

   # plots to check relevance of pluggin estimates against overall average
   hist(obs,type="l",col="steelblue2",nclass=37,proba=T,ylab="",xlab="",main="")
   y=0*xranJ
   r=y
   s=y
   for (j in 1:k){
        y=y+avep[j]*exp(-(xranJ-avmu[j])^2/(2*aves[j]))/sqrt(aves[j])
        r=r+pebest[j]*exp(-(xranJ-mubest[j])^2/(2*sibest[j]))/sqrt(sibest[j])
        s=s+pebesM[j]*exp(-(xranJ-mubesM[j])^2/(2*sibesM[j]))/sqrt(sibesM[j])
	}
   y=y/sqrt(2*pi)
   r=r/sqrt(2*pi)
   s=s/sqrt(2*pi)
   lines(xranJ,y,col="tomato1",lwd=3.7,lty=3)
   lines(xranJ,r,col="chocolate",lwd=3.2,lty=2)
   lines(xranJ,s,col="forestgreen",lwd=2.8,lty=4)
   lines(xranJ,plodens,col="gold",lwd=2)
   legend(.8*max(xranJ),1.5*max(y,r,s),legend=c("mean","MAP","MMAP","ave. dens."),
	lty=c(3,2,4,2),col=c("tomato1","chocola te","forestgreen","gold"),lwd=c(3.7,3.2,2.8,2))

  # saving in file DICxxx.k.ps
   dev.copy(postscript,
   file=paste("DIC",format(runif(1),digits=4),".",k,".ps",sep=""));dev.off()

   # DIC1 (which should not be implemented because of its nonsensical behaviour)
   matlik=apply((avep%*%uno)*exp(-0.5*(dos%*%t(obs)-avmu%*%uno)^2/aves)/
	sqrt(2*pi*aves%*%uno),2,sum)
   pDi1=DICobs+2*sum(log(matlik))
   Di1=DICobs+pDi1

   ## THIS IS DIC2 in the typology 
   # pluggin = Marg MAP
   matlik=apply((pebesM%*%uno)*exp(-0.5*(dos%*%t(obs)-mubesM%*%uno)^2/sibesM)/
   	sqrt(2*pi*sibesM%*%uno),2,sum) 
   pDma=DICobs+2*sum(log(matlik))
   DIC2=DICobs+pDma

   ## THIS IS DIC3 in the typology 
   # w/o pluggin estimate which 
   # is replaced by the average of the current densities
   pDe=DICobs+2*sum(log(avdens))
   DICemp=DICobs+pDe

   ## THIS IS DIC4 in the typology 
   pDt=DICAt-DICB
   DICcompt=DICAt+pDt

   ## THIS IS DIC5 in the typology
   logcomp=2*log(pebest[zbest])-(obs-mubest[zbest])^2/(sibest[zbest])
   	-log(sibest[zbest])
   pD5=DICAt+sum(logcomp)-2*demlog2pi
   DIC5=DICAt+pD5

   ## THIS IS DIC6 in the typology
   # with EM links
   matlik=apply((pebesM%*%uno)*exp(-0.5*(dos%*%t(obs)-mubesM%*%uno)^2/sibesM)/
   	sqrt(2*pi*sibesM%*%uno),2,sum) 
   numera=apply(xlogx((pebesM%*%uno)*exp(-0.5*(dos%*%t(obs)-mubesM%*%uno)^2/
   		(sibesM%*%uno))/sqrt((2*pi)*(sibesM%*%uno))),2,sum)
   pDem=DICAt+2*sum(numera/matlik)
   DICem=DICAt+pDem

   ## THIS IS DIC7 in the typology 
   # DIC7 : use the overall MAP estimator
   comp=exp(-(obs-mubest[zbest])^2/(2*sibest[zbest]))/sqrt(sibest[zbest])
   comp=comp/sqrt(2*pi)
   pD7=DIC1+2*sum(log(comp))
   DIC7=DIC1+pD7

   # Osolete________________________________________________
   # DIC=4 : most likely Z's and estimated theta's
   # The "jitter" avoids ties
   ## THIS IS also DIC6 in the typology 
   #avmz=apply(jitter(avmz,amount=.01),2,rank)
   #zego=0*obs
   #for (j in 1:k)
   #	zego=zego+j*(avmz[j,]==k)
   #logcomp=(obs-avmu[zego])^2/aves[zego]+log(aves[zego])
   #pD4=DIC1-sum(logcomp)-2*demlog2pi
   #DIC4=DIC1+pD4

   ## THIS IS DIC8 in the typology
   pD8=DIC1-DICE
   DIC8=DIC1+pD8
   
   # final preparations
   pdes=c(pDi1,pDma,pDe,pDt,pD5,pDem,pD7,pD8)
   dices=c(Di1,DIC2,DICemp,DICcompt,DIC5,DICem,DIC7,DIC8)
   emthe=c(emresu$p,emresu$mu,emresu$sig)

   list(meth=c("obs|mean","obs|MMAP","obs|emp",
   "comp|RB","comp|MMAP","comp|EM","cond|MAP","cond|condE"),
   typo=c("DIC1","DIC2","DIC3","DIC4","DIC5","DIC6","DIC7","DIC8"),
   pd=pdes,DIC=dices,bestz=zbest,beste=c(pebest,mubest,sibest),ave=c(avep,avmu,aves),
   bestM=c(pebesM,mubesM,sibesM),mle=emthe,beslikes=c(best,besM,compas))
}

lprior=function(p,mu,sig,alpha,esse2,enne,nu,xi){
# 2*log prior
2*sum((alpha-1)*log(p))-
sum(enne*(mu-xi)^2/sig)-
sum(esse2/sig)-sum((nu+1)*log(sig))
}

maktabl=function(litz,waltz,k){
# litz is for list
# waltz is for what
  for (t in 1:length(waltz)){
    j=waltz[t]
    if (j==3)
      digs=format(litz$pd,digits=3)
    if (j==4)
      digs=format(litz$DIC,digits=3)
    if (j==6)
      digs=format(litz$beste,digits=3)
    if (j==7)
      digs=format(litz$ave,digits=3)
    if (j==8)
      digs=format(litz$bestM,digits=3)
    if (j<5){
      dews=1:7
      for (i in 1:6) dews[i]=paste("&",digs[i],sep="")
        dews[7]=paste("&",digs[7],"\\",sep="")
    }
    if (j>4){
      dews=1:(3*k)
      for (i in 1:(3*k-1)) dews[i]=paste("&",digs[i],sep="")
        dews[3*k]=paste("&",digs[3*k],"\\",sep="")
    }
    print(dews,quote=FALSE)
    }
}

xlogx=function(x){
 x*log(x+1e-80)
 }

EMalgo=function(sampl,k,nRep){
# EM algorithm for the k component mixture model
# nRep = number of repetitions of EM (to kill the
# starting value effect)

  nobs=length(sampl)
  uno=t(rep(1,nobs))
  dos=rep(1,k)

  hiya=-10000
  bestmu=rep(0,k)
  bestsi=rep(1,k)
  bestpe=rep(1,k)/k

  for (i in 1:nRep){

    oldlike=-2000000

    # Random initialisation :
    gu=mean(sampl)+1.5*sd(sampl)*rnorm(k)
    su=-var(sampl)/log(runif(k))
    pu=runif(k)
    pu=pu/sum(pu)

    for (t in 2:100){ # 100 is the number of EM steps

    # Allocation probabilities
       probs=(pu%*%uno)*exp(-0.5*(dos%*%t(sampl)-gu%*%uno)^2/su)/sqrt(su%*%uno)
       zeds=t(t(probs)/apply(probs,2,sum)) 

    # EM step
       pu=apply(zeds,1,sum)/nobs
       gu=(sampl%*%t(zeds)/apply(zeds,1,sum))[1,]
       su=apply(zeds*(dos%*%t(sampl)-gu%*%uno)^2,1,sum)/apply(zeds,1,sum)

    # Log-likelihood 
      like=sum(log(apply(probs,2,sum)))

    # Stopping rule
      if (oldlike==like)
        break()
      oldlike=like
      }

    # Update
    if (like>hiya){

      bestmu=gu
      bestsi=su
      bestpe=pu
      hiya=like
      }
  }

  list(mu=gu,sig=su,prob=pu,like=hiya)
}