# A bare-bones S-Plus program for the mixture of exponentials 

# My remission time data file data.dat looks like this (without the # signs): 

#  9.03  34.65  21.47  25.30  80.91  
#
#  9.34   3.88  13.54  15.29  12.10  3.74  23.88   23.87   5.65   1.99   
#
# 27.83  18.44  32.61  19.76  26.68  14.31  21.41   2.01  28.59  10.94
#  6.62   6.41  10.65  62.43  18.50  51.38   1.83   9.24  13.68  23.84
#  3.61   6.45 114.37   8.80   2.55   4.19   2.15  29.15 228.32   5.10
#  0.61  11.73  45.73   8.47  13.43  12.86  12.68  15.65  34.45  19.85
# 46.78  21.94  17.67   1.14   3.08   0.62  27.25  11.19  68.66   2.16


# input data
 x_scan("data.dat")
 r1_sum(x[1:5]); r0_sum(x[6:15]); x_sort(x[16:65]); n_50

# MC sample size and burn-in, as usual
 nmc_5000;  nit_500

# start Gibbs with initial classification and values
# plus vectors to store sampled values 
  mu0.s_mu1.s_NULL
  prob_.5; prob.s_NULL
  z_c(rep(0,25),rep(1,25))
  mu0_mean(x[1:25])+r0/5; mu1_max( mu0+0.1, mean(x[26:50])+r1/10 ) 
  f0_1/mu0; f1_1/mu1
  probistar_rep(0,50)
 
# Gibbs iterations ...
     for (it in 1:(nmc+nit)) {
       # simulate mu0=1/f0 where f0 is Gamma truncated to f0>f1
       t0_sum(1-z); s0_sum((1-z)*x) 
       a_9+t0; b_r0+s0
          u_1-pgamma(b*f1,a); u_min(max(1e-6,u),1-1e-6)
          f0_qgamma(1-runif(1)*u,a)/b; mu0_1/f0



       # simulate mu1=1/f1 where f1 is Gamma truncated to f1<f0
       t1_sum(z); s1_sum(z*x) 
       a_4+t1; b_r1+s1
          u_1-pgamma(b*f0,a); u_min(max(1e-6,u),1-1e-6)
          u_u+runif(1)*(1-u)
          f1_qgamma(1-u,a)/b; mu1_1/f1

       mu1.s_c(mu1.s,mu1)
       mu0.s_c(mu0.s,mu0)

       # simulate (pi) prob given z
       prob_rbeta(1,t1+1,t0+1); prob.s_c(prob.s,prob)

       # simulate new classification indicators z ....
       odds _ ((1-prob)/prob)*(mu1/mu0)*exp(x*(f1-f0))
       pz_1/(1+odds); z_runif(50)<=pz
       if (it > nit) probistar_probistar+ pz
	
#	print(c(prob,mu0,mu1))     
 }

 mu0.s_mu0.s[-(1:nit)]; mu1.s_mu1.s[-(1:nit)]; prob.s_prob.s[-(1:nit)]

 par(mfrow=c(3,3),bty='n')

# posterior means of classication probs ..
  plot(x,probistar/nmc,xlab="ordered data",ylim=c(0,1),
        ylab="posterior classification probs")

# usual plots ... convergence checks
   tsplot(mu0.s,xlab="MCMC iteration",ylab="mu0")
   tsplot(mu1.s,xlab="MCMC iteration",ylab="mu1")
   tsplot(prob.s,xlab="MCMC iteration",ylab="pi")

# margins ...
   hist(mu0.s,nclass=17,probability=T,xlab='mu0')
   hist(mu1.s,nclass=17,probability=T,xlab='mu1')
   hist(prob.s,nclass=17,probability=T,xlab='pi')