"bcdat.rad"<-
structure(.Data = c(45, 6, 0, 46, 46, 7, 17, 7, 37, 0, 4, 15, 11, 22, 46, 46, 
	25, 46, 26, 46, 27, 36, 46, 36, 37, 40, 17, 46, 11, 38, 5, 37, 0, 18, 
	24, 36, 5, 19, 17, 24, 32, 33, 19, 37, 34, 36, Inf, 10, 7, Inf, Inf, 16,
	Inf, 14, 44, 8, 11, Inf, 15, Inf, Inf, Inf, 37, Inf, 40, Inf, 34, 44, 
	Inf, 48, Inf, Inf, 25, Inf, 18, Inf, 12, Inf, 5, Inf, Inf, Inf, 11, 35, 
	25, Inf, Inf, Inf, 26, Inf, Inf, Inf), .Dim = c(46, 2), .Dimnames = 
	list(c("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"), c("left", 
	"right")))
"bcdat.radchem"<-
structure(.Data = c(8, 0, 24, 17, 17, 24, 16, 13, 11, 16, 18, 17, 32, 23, 44, 
	14, 0, 5, 12, 11, 33, 31, 13, 19, 34, 13, 16, 35, 15, 11, 22, 10, 30, 
	13, 10, 8, 4, 11, 14, 4, 34, 30, 18, 16, 35, 21, 11, 48, 12, 22, 31, 27,
	23, 30, 24, Inf, 13, 20, 25, 26, Inf, Inf, 48, 17, 5, 8, 20, Inf, 40, 
	Inf, 39, 32, Inf, Inf, 24, Inf, 22, 17, 32, 35, 34, Inf, 17, 21, 9, Inf,
	19, 8, Inf, 36, 24, 60, 39, Inf, 20, Inf), .Dim = c(48, 2), .Dimnames
	 = list(c("47", "48", "49", "50", "51", "52", "53", "54", "55", "56", 
	"57", "58", "59", "60", "61", "62", "63", "64", "65", "66", "67", "68", 
	"69", "70", "71", "72", "73", "74", "75", "76", "77", "78", "79", "80", 
	"81", "82", "83", "84", "85", "86", "87", "88", "89", "90", "91", "92", 
	"93", "94"), c("left", "right")))
"cisnrxx"<-
function(L, ci, func, nrep, npass = 3, min.size = 500, tol = 0, reuse = F, seed
	 = reuse)
{
# L is output from gibbs1.
	if(is.null(nrep)) nrep <- pmax(min.size, ceiling(L$nrep/2^((npass - 1):
			0)))
	else if(length(nrep)==1) nrep<-rep(nrep,npass)
	else npass<-length(nrep)
	if(length(seed) == 3)
		seed <- seed
	else if(seed)
		seed <- ceiling(ranf(3) * c(32, 9999, 9999))
	else seed <- NULL
	start.ranf(seed)
	tlim <- intvlxx(L, func, level = ci, nrep = nrep[1], new = reuse, tol
		 = tol)
	if(npass > 1) for(i in 2:npass) {
			if(reuse)
				start.ranf(seed = seed)
			tlim <- intvlxx(L, func, level = ci, nrep = nrep[i], 
				guess = tlim$limits, new = reuse, tol = tol)
		}
# Convert upper limits to proper values.
	tlim$limits[2,  ] <- 1 - tlim$limits[2,  ]
	c(tlim,list(nrep.ci=nrep))
}
"cmpall"<-
function(L, domcse = F)
{
# This function compares the survival functions and
# standard errors from different gibbs1, sis1, sis2  or ritcen
# runs. The input is a list each of whose elements is an output
# list from gibbs1, sis1, sis2 or ritcen.
	cmp <- function(s1, s2, nam)
	{
		comp <- function(x, y)
		{
			max(abs(x - y))
		}
		n <- length(nam)
		ans <- rep(NA, n)
		timx <- s1$time
		timy <- s2$time
		timx <- timx[ - length(timx)]
		timy <- timy[ - length(timy)]
		sub <- intersect(timx, timy)
		if(length(sub) == 0)
			return(ans)
		subx <- match(sub, timx)
		suby <- match(sub, timy)
		for(i in 1:n) {
			x <- s1[[nam[i]]]
			y <- s2[[nam[i]]]
			if((length(x) <= 1) || (length(y) <= 1))
				ans[i] <- NA
			else ans[i] <- comp(x[subx], y[suby])
		}
		ans
	}
	nam <- c("surv", "sigma")
	if(domcse)
		nam <- c(nam, "semu")
	n <- length(L)
	ans <- c()
	for(i in 1:(n - 1))
		for(j in (i + 1):n) {
			ans <- rbind(ans, c(i, j, cmp(L[[i]], L[[j]], nam)))
		}
	dimnames(ans) <- list(NULL, c("fit1", "fit2", nam))
	ans
}
"gbs1ci.call"<-
function(L)
{
# The components of the list L must be ordered as in the .Fortran
# function call below.
	storage.mode(L[[1]]) <- "double"
	storage.mode(L[[2]]) <- "double"
	storage.mode(L[[3]]) <- "integer"
	storage.mode(L[[4]]) <- "integer"
	storage.mode(L[[5]]) <- "integer"
	storage.mode(L[[6]]) <- "double"
	storage.mode(L[[7]]) <- "double"
	storage.mode(L[[8]]) <- "double"
	storage.mode(L[[9]]) <- "integer"
	storage.mode(L[[10]]) <- "integer"
	storage.mode(L[[11]]) <- "integer"
	storage.mode(L[[12]]) <- "integer"
	storage.mode(L[[13]]) <- "double"
	storage.mode(L[[14]]) <- "integer"
	storage.mode(L[[15]]) <- "integer"
	storage.mode(L[[16]]) <- "double"
	storage.mode(L[[17]]) <- "integer"
	storage.mode(L[[18]]) <- "logical"
	storage.mode(L[[19]]) <- "double"
	storage.mode(L[[20]]) <- "double"
	storage.mode(L[[21]]) <- "double"
	.Fortran("gbs1ci",
		lo = L[[1]],
		hi = L[[2]],
		inlo = L[[3]],
		inhi = L[[4]],
		nn = L[[5]],
		aa0 = L[[6]],
		bb0 = L[[7]],
		totm0 = L[[8]],
		nrep = L[[9]],
		g0 = L[[10]],
		g1 = L[[11]],
		gextra = L[[12]],
		cc = L[[13]],
		intim = L[[14]],
		nc = L[[15]],
		time = L[[16]],
		ntim = L[[17]],
		new = L[[18]],
		tol = L[[19]],
		guess = L[[20]],
		beta = L[[21]])
}
"gibbs1"<-
function(data, prior, nrep, control = c(10, 1, 3), seed = NULL, eps = 0.01, 
	adtim = NULL, domcse = T, ci = 0, nbatch = 20, new = T, npass.ci = 3, 
	nrep.ci = NULL, tol = 0, reuse = F)
{
#
# Required arguments:
# ------------------
# data is a matrix containing the data intervals.  The matrix has two
#   columns.  Each row of the matrix gives one data interval.  The time
#   value infinity may be coded as NA, Inf, or -99.
# prior is a vector specifying the prior distribution.  It contains
#   three values a, b, totm0.  The values a, b specify the gamma mixing
#   distribution.  The value totm0 is the total mass alpha(R) of the 
#   Dirichlet measure.
# nrep is the number of iterations of the Gibbs sampler to record.
#
# Optional arguments:
# ------------------
# control is a vector of three values (go, g1, gextra) which control
#    the running of the Gibbs sampler.
#    g0 is the warm-up or burn-in period.  The first g0-1 iterations of
#      the "complete" Gibbs sampler are discarded.
#    g1 is the spacing.  Every g1 iterations of the "complete" Gibbs
#      sampler, we record various statistics (keep running totals).
#    gextra controls the spacing of the "extra step" which randomizes
#      the position of the atoms.  The extra step is inserted every
#      gextra sweeps of the "basic" Gibbs sampler.  The gextra sweeps
#      of the "basic" Gibbs sampler plus the extra step make up one
#      iteration of the "complete" Gibbs sampler.
# seed is a vector of three integer seeds for the random number 
#    generators.  The first seed should be between 1 and 32.
# eps is a small positive value.  The program currently does not
#    handle exact observed death times t, but replaces these by 
#    interval censored observations (t-eps,t).
# adtim is a vector of additional time points at which the survival
#    function is to be computed.  If this is not specified, then
#    the survival function is only computed at the endpoints of the
#    data intervals.
#
	infty <- -99
	tmp <- lohi(data, eps)
	lo <- tmp$lo
	hi <- tmp$hi
	nn <- length(lo)
	aa0 <- prior[1]
	bb0 <- prior[2]
	totm0 <- prior[3]
	g0 <- control[1]
	g1 <- control[2]
	gextra <- control[3]
	start.ranf(seed)
	misc <- setup(lo, hi, adtim)
	inlo <- misc$inlo
	inhi <- misc$inhi
	cc <- misc$cc
	nc <- length(cc)
	intim <- misc$intim
	time <- misc$time
	ntim <- misc$ntim
	surv <- rep(0, ntim)
	sigma <- rep(0, ntim)
	muth <- 0
	varth <- 0
	L <- list(lo = lo, hi = hi, inlo = inlo, inhi = inhi, nn = nn, aa0 = 
		aa0, bb0 = bb0, totm0 = totm0, nrep = nrep, g0 = g0, g1 = g1, 
		gextra = gextra, cc = cc, intim = intim, nc = nc, time = time, 
		ntim = ntim, new = new, surv = surv, sigma = sigma, muth = muth,
		varth = varth, doci = F, tol = tol, guess = 0, beta = 0)
	if(domcse & (nbatch > 1)) {
		ntot <- nbatch * (L$nrep <- ceiling(L$nrep/nbatch))
		ans <- gibbs1.call(L)
		L$new <- F
		L$g0 <- if(length(control) == 4) control[4] else L$g1
		muth <- ans$mu
		varth <- ans$varth
		semu <- (av.surv <- ans$surv)^2
		sz2 <- (sz <- (s.sigma <- ans$sigma) - av.surv^2)^2
		for(i in 2:nbatch) {
			ans <- gibbs1.call(L)
			muth <- muth + ans$muth
			varth <- varth + ans$varth
			srv <- ans$surv
			sig <- ans$sigma
			av.surv <- av.surv + srv
			semu <- semu + srv^2
			s.sigma <- s.sigma + sig
			sz <- sz + (z <- sig - srv^2)
			sz2 <- sz2 + z^2
		}
		ans$muth <- muth/nbatch
		ans$varth <- varth/nbatch
		ans$surv <- (av.surv <- av.surv/nbatch)
		semu <- (pmax((semu/nbatch - av.surv^2)/(nbatch - 1), 0))^0.5
		ans$sigma <- (pmax(s.sigma/nbatch - (av.surv)^2, 0))^0.5
		sesig <- ((sz2/nbatch - (sz/nbatch)^2)/(nbatch - 1))^0.5/(2 * 
			ans$sigma)
		if(time[1] == 0)
			sesig[1] <- 0
		ans$nrep <- ntot
		ans$g0 <- g0
		ans <- c(ans, list(semu = semu, sesig = sesig, nbatch = nbatch)
			)
	}
	else {
		ans <- gibbs1.call(L)
		ans$sigma <- (pmax(ans$sigma - (ans$surv)^2, 0))^0.5
	}
	ans$varth <- ans$varth - (ans$muth)^2
	if(ci > 0)
		ans <- c(ans, ci = ci, cisnrxx(ans[1:22], ci, "gibbs1.call",
			nrep.ci, npass.ci, min.size = 500, tol = tol, reuse = 
			reuse))
	class(ans) <- "mixdir"
	ans
}
"gibbs1.call"<-
function(L)
{
# The components of the list L must be ordered as in the .Fortran
# function call below.
	storage.mode(L[[1]]) <- "double"
	storage.mode(L[[2]]) <- "double"
	storage.mode(L[[3]]) <- "integer"
	storage.mode(L[[4]]) <- "integer"
	storage.mode(L[[5]]) <- "integer"
	storage.mode(L[[6]]) <- "double"
	storage.mode(L[[7]]) <- "double"
	storage.mode(L[[8]]) <- "double"
	storage.mode(L[[9]]) <- "integer"
	storage.mode(L[[10]]) <- "integer"
	storage.mode(L[[11]]) <- "integer"
	storage.mode(L[[12]]) <- "integer"
	storage.mode(L[[13]]) <- "double"
	storage.mode(L[[14]]) <- "integer"
	storage.mode(L[[15]]) <- "integer"
	storage.mode(L[[16]]) <- "double"
	storage.mode(L[[17]]) <- "integer"
	storage.mode(L[[18]]) <- "logical"
	storage.mode(L[[19]]) <- "double"
	storage.mode(L[[20]]) <- "double"
	storage.mode(L[[21]]) <- "double"
	storage.mode(L[[22]]) <- "double"
	storage.mode(L[[23]]) <- "logical"
	storage.mode(L[[24]]) <- "double"
	storage.mode(L[[25]]) <- "double"
	storage.mode(L[[26]]) <- "double"
	.Fortran("gibbs1",
		lo = L[[1]],
		hi = L[[2]],
		inlo = L[[3]],
		inhi = L[[4]],
		nn = L[[5]],
		aa0 = L[[6]],
		bb0 = L[[7]],
		totm0 = L[[8]],
		nrep = L[[9]],
		g0 = L[[10]],
		g1 = L[[11]],
		gextra = L[[12]],
		cc = L[[13]],
		intim = L[[14]],
		nc = L[[15]],
		time = L[[16]],
		ntim = L[[17]],
		new = L[[18]],
		surv = L[[19]],
		sigma = L[[20]],
		muth = L[[21]],
		varth = L[[22]],
		doci = L[[23]],
		tol = L[[24]],
		guess = L[[25]],
		beta = L[[26]])
}
"gibbs2"<-
function(data, prior, nrep, control = c(10, 1, 3), seed = NULL, eps = 0.01, 
	adtim = NULL, domcse = T, nbatch = 20, new = T, ord = T)
{
#
# Required arguments:
# ------------------
# data is a matrix containing the data intervals.  The matrix has two
#   columns.  Each row of the matrix gives one data interval.  The time
#   value infinity may be coded as NA, Inf, or -99.
# prior is a vector specifying the prior distribution.  It contains
#   three values a, b, totm0.  The values a, b specify the gamma mixing
#   distribution.  The value totm0 is the total mass alpha(R) of the 
#   Dirichlet measure.
# nrep is the number of iterations of the Gibbs sampler to record.
#
# Optional arguments:
# ------------------
# control is a vector of three values (go, g1, gextra) which control
#    the running of the Gibbs sampler.
#    g0 is the warm-up or burn-in period.  The first g0-1 iterations of
#      the "complete" Gibbs sampler are discarded.
#    g1 is the spacing.  Every g1 iterations of the "complete" Gibbs
#      sampler, we record various statistics (keep running totals).
#    gextra controls the spacing of the "extra step" which randomizes
#      the position of the atoms.  The extra step is inserted every
#      gextra sweeps of the "basic" Gibbs sampler.  The gextra sweeps
#      of the "basic" Gibbs sampler plus the extra step make up one
#      iteration of the "complete" Gibbs sampler.
# seed is a vector of three integer seeds for the random number 
#    generators.  The first seed should be between 1 and 32.
# eps is a small positive value.  The program currently does not
#    handle exact observed death times t, but replaces these by 
#    interval censored observations (t-eps,t).
# adtim is a vector of additional time points at which the survival
#    function is to be computed.  If this is not specified, then
#    the survival function is only computed at the endpoints of the
#    data intervals.
#
	infty <- -99
	tmp <- lohi(data, eps)
	lo <- tmp$lo
	hi <- tmp$hi
	nn <- length(lo)
	aa0 <- prior[1]
	bb0 <- prior[2]
	totm0 <- prior[3]
	g0 <- control[1]
	g1 <- control[2]
	gextra <- control[3]
	if(ord) {
		phi <- bb0^aa0/(bb0 + hi)^aa0
		phi[hi == infty] <- 0
		plo <- bb0^aa0/(bb0 + lo)^aa0
		prb <- plo - phi
		prb[tmp$tie] <- 0
		pi <- order(prb)
		lo <- lo[pi]
		hi <- hi[pi]
		grp <- .Fortran("grnest",
			lo = as.double(lo),
			hi = as.double(hi),
			nn = as.integer(nn),
			ngp = as.integer(0))
		lo <- grp$lo
		hi <- grp$hi
	}
	start.ranf(seed = seed)
	misc <- setup(lo, hi, adtim)
	inlo <- misc$inlo
	inhi <- misc$inhi
	cc <- misc$cc
	nc <- length(cc)
	intim <- misc$intim
	time <- misc$time
	ntim <- misc$ntim
	surv <- rep(0, ntim)
	sigma <- rep(0, ntim)
	mu <- 0
	var <- 0
	L <- list(lo = lo, hi = hi, inlo = inlo, inhi = inhi, nn = nn, aa0 = 
		aa0, bb0 = bb0, totm0 = totm0, nrep = nrep, g0 = g0, g1 = g1, 
		gextra = gextra, cc = cc, intim = intim, nc = nc, time = time, 
		ntim = ntim, new = new, surv = surv, sigma = sigma, mu = mu, 
		var = var)
	if(domcse & (nbatch > 1)) {
		ntot <- nbatch * (L$nrep <- ceiling(L$nrep/nbatch))
		ans <- gibbs2.call(L)
		L$new <- F
		L$g0 <- if(length(control) == 4) control[4] else L$g1
		muth <- ans$mu
		varth <- ans$var
		semu <- (av.surv <- ans$surv)^2
		sz2 <- (sz <- (s.sigma <- ans$sigma) - av.surv^2)^2
		for(i in 2:nbatch) {
			ans <- gibbs2.call(L)
			muth <- muth + ans$mu
			varth <- varth + ans$var
			srv <- ans$surv
			sig <- ans$sigma
			av.surv <- av.surv + srv
			semu <- semu + srv^2
			s.sigma <- s.sigma + sig
			sz <- sz + (z <- sig - srv^2)
			sz2 <- sz2 + z^2
		}
		ans$mu <- muth/nbatch
		ans$var <- varth/nbatch
		ans$surv <- (av.surv <- av.surv/nbatch)
		semu <- (pmax((semu/nbatch - av.surv^2)/(nbatch - 1), 0))^0.5
		ans$sigma <- (pmax(s.sigma/nbatch - (av.surv)^2, 0))^0.5
		sesig <- ((sz2/nbatch - (sz/nbatch)^2)/(nbatch - 1))^0.5/(2 * 
			ans$sigma)
		if(time[1] == 0)
			sesig[1] <- 0
		ans$nrep <- ntot
		ans <- c(ans, list(semu = semu, sesig = sesig, nbatch = nbatch)
			)
	}
	else {
		ans <- gibbs2.call(L)
		ans$sigma <- (pmax(ans$sigma - (ans$surv)^2, 0))^0.5
	}
	ans$var <- ans$var - (ans$mu)^2
	class(ans) <- "mixdir"
	ans
}
"gibbs2.call"<-
function(L)
{
# The components of the list L must be ordered as in the .Fortran
# function call below.
	storage.mode(L[[1]]) <- "double"
	storage.mode(L[[2]]) <- "double"
	storage.mode(L[[3]]) <- "integer"
	storage.mode(L[[4]]) <- "integer"
	storage.mode(L[[5]]) <- "integer"
	storage.mode(L[[6]]) <- "double"
	storage.mode(L[[7]]) <- "double"
	storage.mode(L[[8]]) <- "double"
	storage.mode(L[[9]]) <- "integer"
	storage.mode(L[[10]]) <- "integer"
	storage.mode(L[[11]]) <- "integer"
	storage.mode(L[[12]]) <- "integer"
	storage.mode(L[[13]]) <- "double"
	storage.mode(L[[14]]) <- "integer"
	storage.mode(L[[15]]) <- "integer"
	storage.mode(L[[16]]) <- "double"
	storage.mode(L[[17]]) <- "integer"
	storage.mode(L[[18]]) <- "logical"
	storage.mode(L[[19]]) <- "double"
	storage.mode(L[[20]]) <- "double"
	storage.mode(L[[21]]) <- "double"
	storage.mode(L[[22]]) <- "double"
	.Fortran("gibbs2",
		lo = L[[1]],
		hi = L[[2]],
		inlo = L[[3]],
		inhi = L[[4]],
		nn = L[[5]],
		aa0 = L[[6]],
		bb0 = L[[7]],
		totm0 = L[[8]],
		nrep = L[[9]],
		g0 = L[[10]],
		g1 = L[[11]],
		gextra = L[[12]],
		cc = L[[13]],
		intim = L[[14]],
		nc = L[[15]],
		time = L[[16]],
		ntim = L[[17]],
		new = L[[18]],
		surv = L[[19]],
		sigma = L[[20]],
		mu = L[[21]],
		var = L[[22]])
}
"guess.ci"<-
function(m1, sigma, ci, zero)
{
	if(zero) {
		m1 <- m1[-1]
		sigma <- sigma[-1]
	}
	m2 <- (var <- sigma^2) + m1^2
	a <- (m1 * (m1 - m2))/var
	b <- ((1/m1) - 1) * a
	p <- (1 - ci)/2
	ans <- t(cbind(qbeta(p, a, b), qbeta(p, b, a)))
	if(zero)
		ans <- cbind(c(1, 0), ans)
	ans
}
"improve.guess"<-
function(y, p, tol = 0)
{
# y[1,] is the guess for the confidence limit.
# y[2,] is the estimated lower tail probability at y[1,].
# y[3,] is the estimated density at y[1,].
# y[4,] is the estimated derivative at y[1,].
	nc <- ncol(y)
	.Fortran("impgss",
		as.double(y),
		as.integer(nc),
		as.double(p),
		as.double(tol),
		double(nc))[[5]]
}
"intvlxx"<-
function(L, func, level = 0.94999999999999996, new = F, nrep = L$nrep, guess = 
	NULL, tol = 0)
{
	if(is.null(guess))
		guess <- guess.ci(L$surv, L$sigma, level, L$time[1] == 0)
	storage.mode(guess) <- "single"
	beta <- matrix(0, 6, L$ntim)
	if(func == "gibbs1.call")
		L$new <- new
	L$nrep <- nrep
	beta <- eval(call(func, c(L, list(doci = T, tol = tol, guess = guess, 
		beta = beta))))$beta
	lim <- list(limits = limxx(guess, beta, level, tol = tol))
	c(lim, list(prev.lim = guess, prev.prob = beta[c(1, 4),  ]))
}
"limxx"<-
function(y0, f, level, tol = 0)
{
	p <- (1 - level)/2
	lower <- improve.guess(rbind(y0[1,  ], f[1:3,  ]), p, tol = tol)
	upper <- improve.guess(rbind(y0[2,  ], f[4:6,  ]), p, tol = tol)
	rbind(lower, upper)
}
"lohi"<-
function(data, eps)
{
# This function is used within sis1 and gibbs1.
# Users need not worry about this function.
	infty <- -99	
	# If column 2 of data consists entirely of 1's and 0's  we assume
#    these are indicators of death and right-censoring respectively.
	if(all((data[, 2] == 0) | (data[, 2] == 1))) {
		dead <- data[, 2] == 1
		data[dead, 2] <- data[dead, 1]
		data[!dead, 2] <- infty
	}
	lo <- data[, 1]
	hi <- data[, 2]
	hi[hi == Inf] <- infty
	hi[is.na(hi)] <- infty
	lo[lo ==  - Inf] <- 0
	lo[is.na(lo)] <- 0
	tie <- (lo == hi)
	lo[tie] <- lo[tie] - eps
	list(lo = lo, hi = hi, tie = tie)
}
"parmod"<-
function(data, prior, nrep, control = c(10, 1, 3), seed = NULL, eps = 0.01, 
	adtim = NULL, domcse = T, nbatch = 20, new = T)
{
#
# Required arguments:
# ------------------
# data is a matrix containing the data intervals.  The matrix has two
#   columns.  Each row of the matrix gives one data interval.  The time
#   value infinity may be coded as NA, Inf, or -99.
# prior is a vector specifying the prior distribution.  It contains
#   three values a, b, totm0.  The values a, b specify the gamma mixing
#   distribution.  The value totm0 is the total mass alpha(R) of the 
#   Dirichlet measure.
# nrep is the number of iterations of the Gibbs sampler to record.
#
# Optional arguments:
# ------------------
# control is a vector of three values (go, g1, gextra) which control
#    the running of the Gibbs sampler.
#    g0 is the warm-up or burn-in period.  The first g0-1 iterations of
#      the "complete" Gibbs sampler are discarded.
#    g1 is the spacing.  Every g1 iterations of the "complete" Gibbs
#      sampler, we record various statistics (keep running totals).
#    gextra controls the spacing of the "extra step" which randomizes
#      the position of the atoms.  The extra step is inserted every
#      gextra sweeps of the "basic" Gibbs sampler.  The gextra sweeps
#      of the "basic" Gibbs sampler plus the extra step make up one
#      iteration of the "complete" Gibbs sampler.
# seed is a vector of three integer seeds for the random number 
#    generators.  The first seed should be between 1 and 32.
# eps is a small positive value.  The program currently does not
#    handle exact observed death times t, but replaces these by 
#    interval censored observations (t-eps,t).
# adtim is a vector of additional time points at which the survival
#    function is to be computed.  If this is not specified, then
#    the survival function is only computed at the endpoints of the
#    data intervals.
#
	infty <- -99
	tmp <- lohi(data, 0)
	lo <- tmp$lo
	hi <- tmp$hi
	nn <- length(lo)	
	# Absorb complete observations and right-censored observations
# into the prior.
	aa0 <- prior[1] + sum(tie <- lo == hi)
	bb0 <- prior[2] + sum(lo[rcen <- hi == infty]) + sum(lo[tie])
	del <- tie | rcen
	adtim <- unique(sort(c(adtim, lo[del])))
	lo <- lo[!del]
	hi <- hi[!del]
	if(length(lo) == 0)
		stop("There must be at least one interval censored observation.\n"
			)
	nn <- length(lo)
	g0 <- control[1]
	g1 <- control[2]
	gextra <- control[3]
	start.ranf(seed = seed)
	misc <- setup(lo, hi, adtim)
	time <- misc$time
	ntim <- misc$ntim
	surv <- rep(0, ntim)
	sigma <- rep(0, ntim)
	mu <- 0
	var <- 0
	L <- list(lo = lo, hi = hi, nn = nn, aa0 = aa0, bb0 = bb0, nrep = nrep, 
		g0 = g0, g1 = g1, gextra = gextra, time = time, ntim = ntim, 
		new = new, surv = surv, sigma = sigma, mu = mu, var = var)
	if(domcse & (nbatch > 1)) {
		ntot <- nbatch * (L$nrep <- ceiling(L$nrep/nbatch))
		ans <- parmod.call(L)
		L$new <- F
		L$g0 <- if(length(control) == 4) control[4] else L$g1
		muth <- ans$mu
		varth <- ans$var
		semu <- (av.surv <- ans$surv)^2
		sz2 <- (sz <- (s.sigma <- ans$sigma) - av.surv^2)^2
		for(i in 2:nbatch) {
			ans <- parmod.call(L)
			muth <- muth + ans$mu
			varth <- varth + ans$var
			srv <- ans$surv
			sig <- ans$sigma
			av.surv <- av.surv + srv
			semu <- semu + srv^2
			s.sigma <- s.sigma + sig
			sz <- sz + (z <- sig - srv^2)
			sz2 <- sz2 + z^2
		}
		ans$mu <- muth/nbatch
		ans$var <- varth/nbatch
		ans$surv <- (av.surv <- av.surv/nbatch)
		semu <- (pmax((semu/nbatch - av.surv^2)/(nbatch - 1), 0))^0.5
		ans$sigma <- (pmax(s.sigma/nbatch - (av.surv)^2, 0))^0.5
		sesig <- ((sz2/nbatch - (sz/nbatch)^2)/(nbatch - 1))^0.5/(2 * 
			ans$sigma)
		if(time[1] == 0)
			sesig[1] <- 0
		ans$nrep <- ntot
		ans$g0 <- g0
		ans <- c(ans, list(semu = semu, sesig = sesig, nbatch = nbatch)
			)
	}
	else {
		ans <- parmod.call(L)
		ans$sigma <- (pmax(ans$sigma - (ans$surv)^2, 0))^0.5
	}
	ans$var <- ans$var - (ans$mu)^2
	class(ans) <- "mixdir"
	ans
}
"parmod.call"<-
function(L)
{
# The components of the list L must be ordered as in the .Fortran
# function call below.
	storage.mode(L[[1]]) <- "double"
	storage.mode(L[[2]]) <- "double"
	storage.mode(L[[3]]) <- "integer"
	storage.mode(L[[4]]) <- "double"
	storage.mode(L[[5]]) <- "double"
	storage.mode(L[[6]]) <- "integer"
	storage.mode(L[[7]]) <- "integer"
	storage.mode(L[[8]]) <- "integer"
	storage.mode(L[[9]]) <- "integer"
	storage.mode(L[[10]]) <- "double"
	storage.mode(L[[11]]) <- "integer"
	storage.mode(L[[12]]) <- "logical"
	storage.mode(L[[13]]) <- "double"
	storage.mode(L[[14]]) <- "double"
	storage.mode(L[[15]]) <- "double"
	storage.mode(L[[16]]) <- "double"
	.Fortran("parmod",
		lo = L[[1]],
		hi = L[[2]],
		nn = L[[3]],
		aa0 = L[[4]],
		bb0 = L[[5]],
		nrep = L[[6]],
		g0 = L[[7]],
		g1 = L[[8]],
		gextra = L[[9]],
		time = L[[10]],
		ntim = L[[11]],
		new = L[[12]],
		surv = L[[13]],
		sigma = L[[14]],
		mu = L[[15]],
		var = L[[16]])
}
"plot.mixdir"<-
function(L, band = T, mult = 1, xlab = "time", ylab = "surv", new = T, lty1 = 1,
	lty2 = 2, use.ci = T, ...)
{
# Set band = F if you do NOT wish to print the band
#    about the survival function.
# The band is plus or minus mult*sigma on either side
#    of the survival function (where sigma is the posterior
#    standard deviation).  The default is mult = 1, but mult
#    can be set to any value.
# If you want to add to an existing plot, set new=F.
# The values lty1 and lty2 control the type of line plotted for
#    the survival function and bands, respectively.
#----------------------------------------------------
	x <- L$time[1:L$ntim]
	y <- L$surv
	if(band) {
		if(use.ci && !is.null(L$limits)) {
			lower <- (L$limits)[1,  ]
			upper <- (L$limits)[2,  ]
		}
		else {
			se <- L$sigma
			upper <- y + mult * se
			lower <- y - mult * se
		}
		yp <- c(0, 1, y, upper, lower)
		if(new)
			plot(cbind(range(x), range(yp)), type = "n", xlab = 
				xlab, ylab = ylab, ...)
		lines(x, y, lty = lty1)
		lines(x, upper, lty = lty2)
		lines(x, lower, lty = lty2)
	}
	else {
		if(new)
			plot(range(x), c(0, 1), type = "n", xlab = xlab, ylab
				 = ylab, ...)
		lines(x, y, lty = lty1)
	}
}
"ritcen"<-
function(data, prior, nrep, seed = NULL, eps = 0.01, adtim = NULL, domcse = T, 
	nthint = 10, ci = 0, npass.ci = 3, nrep.ci = NULL, tol = 0, reuse = F)
{
#-----------------------------------------------------
# Note: This function can only handle right-censoring.
#
# Required arguments:
# ------------------
# data is a matrix containing the data intervals.  The matrix has two
#   columns.  Each row of the matrix gives one observation or data interval.
#   The data may be in one of two different formats.  
#  
#  Format 1: Each row gives the lower and upper endpoints of the time
#    interval for that observation.  An observed death time (an uncensored)
#    observation) is represented as a degenerate interval with both
#    endpoints equal.  For right-censored observations, the time
#    value infinity may be coded as NA, Inf, or -99.  This is the same
#    format used by gibbs1 and sis1.
#
#  Format 2: Column 1 gives the times, and column 2 gives an indicator of
#    status with 1 = dead (uncensored) and 0 = alive (right-censored). This
#    format is suitable for passing to the Splus survfit function.
#
#  The program decides which format is being used by looking at the second
#  column of data.  If it consists entirely of 0's and 1's, then format 2
#  is assumed.
# 
# prior is a vector specifying the prior distribution.  It contains
#   three values a, b, totm.  The values a, b specify the gamma mixing
#   distribution.  The value totm is the total mass alpha(R) of the 
#   Dirichlet measure.
# nrep is the simulation sample size.  The final answers are averages
#   over nrep i.i.d. repetitions.
#
# Optional arguments:
# ------------------
# seed is a vector of three integer seeds for the random number 
#    generators.  The first seed should be between 1 and 32.
# eps is a small positive value.  For convenience, uncensored 
#    observations (at time t) are viewed as belonging to the
#    interval (t-eps,t).  Computation and plotting of the posterior
#    mean and standard deviation of the survival function will
#    be done at both the values t and t-eps.
# adtim is a vector of additional time points at which the survival
#    function is to be computed.  If this is not specified, then
#    the survival function is only computed at the endpoints of the
#    data intervals.
# domcse is a logical variable (T or F)  specifying whether or not
#    Monte Carlo standard errors for the posterior mean and variance
#    will be computed.
# nthint is the number of theta intervals used in the rejection scheme
#    for generating from the posterior distribution of theta.
#    The average number of variates generated (until one is accepted)
#    is recorded in the component "eff" of the answer list.  If eff
#    gets large (say, close to 2), then nthint should be increased.
#
#
# Returned Value (partial description)
# --------------
# time is a list of the time points at which the quantities below
#    are calculated.
# surv is an estimated of the posterior mean of the survival function
#    at the various time points.
# semu is an estimate of the simulation standard error of surv.  This
#    tells us the accuracy of surv.  If greater accuracy is desired, 
#    you must increase the value of nrep.
# sigma is an estimate of the posterior standard deviation of the survival
#    function.
# sesig is an estimate of the simulation standard error of sigma.
#---------------------------------------------------
	infty <- -99
	errmess <- "This program can only handle right-censoring.\n"
	tie <- -33
	nn <- nrow(data)
	z <- data[, 1]
	if(!(all(is.finite(z)) && all(z >= 0))) {
		cat("Bad values in column 1 of data.\n")
		cat(errmess)
		return(NULL)
	}
	if(all((data[, 2] == 0) | (data[, 2] == 1))) {
		delt <- data[, 2]
	}
	else {
		hi <- data[, 2]
		hi[(hi == infty) | is.na(hi)] <- Inf
		if(any(sub <- (hi == tie)))
			hi[sub] <- z[sub]
		if(any(z > hi) || any(is.finite(hi) & (hi > z))) {
			cat("Bad values in column 2 of data.\n")
			cat(errmess)
			return(NULL)
		}
		delt <- 1 * (z == hi)
	}
	aa <- prior[1]
	bb <- prior[2]
	totm <- prior[3]
	start.ranf(seed = seed)
	pi <- rev(order(z,  - delt))
	z <- z[pi]
	delt <- delt[pi]
	uz <- unique(sort(z))
	eps <- min(eps, 0.5 * min(uz[-1] - uz[ - length(uz)]))
	time <- unique(sort(c(z, z[delt == 1] - eps, adtim)))
	ntim <- length(time)
	intim <- match(z, time) - delt
	time <- c(time, infty)
	sigma <- surv <- rep(0, ntim)
	if(domcse)
		semu <- sesig <- surv
	else semu <- sesig <- 0
	nadtim <- length(adtim)
	muth <- 0
	varth <- 0
	eff <- 0
	L <- list(z = z, delt = delt, intim = intim, nn = nn, time = time, ntim
		 = ntim, aa = aa, bb = bb, totm = totm, nrep = nrep, domcse = 
		domcse, nthint = nthint, surv = surv, sigma = sigma, semu = 
		semu, sesig = sesig, muth = muth, varth = varth, eff = eff, 
		doci = F, tol = tol, guess = 0, beta = 0)
	ans <- ritcen.call(L)
	if(ci > 0)
		ans <- c(ans, ci = ci, cisnrxx(ans[1:19], ci, "ritcen.call",
			nrep.ci, npass.ci, min.size = 500, tol = tol, reuse = 
			reuse))
	if(time[1] == 0)
		ans$sesig[1] <- 0
	class(ans) <- "mixdir"
	ans
}
"ritcen.call"<-
function(L)
{
	storage.mode(L[[1]]) <- "double"
	storage.mode(L[[2]]) <- "integer"
	storage.mode(L[[3]]) <- "integer"
	storage.mode(L[[4]]) <- "integer"
	storage.mode(L[[5]]) <- "double"
	storage.mode(L[[6]]) <- "integer"
	storage.mode(L[[7]]) <- "double"
	storage.mode(L[[8]]) <- "double"
	storage.mode(L[[9]]) <- "double"
	storage.mode(L[[10]]) <- "integer"
	storage.mode(L[[11]]) <- "logical"
	storage.mode(L[[12]]) <- "integer"
	storage.mode(L[[13]]) <- "double"
	storage.mode(L[[14]]) <- "double"
	storage.mode(L[[15]]) <- "double"
	storage.mode(L[[16]]) <- "double"
	storage.mode(L[[17]]) <- "double"
	storage.mode(L[[18]]) <- "double"
	storage.mode(L[[19]]) <- "double"
	storage.mode(L[[20]]) <- "logical"
	storage.mode(L[[21]]) <- "double"
	storage.mode(L[[22]]) <- "double"
	storage.mode(L[[23]]) <- "double"
	.Fortran("ritcen",
		z = L[[1]],
		delt = L[[2]],
		intim = L[[3]],
		nn = L[[4]],
		time = L[[5]],
		ntim = L[[6]],
		aa = L[[7]],
		bb = L[[8]],
		totm = L[[9]],
		nrep = L[[10]],
		domcse = L[[11]],
		nthint = L[[12]],
		surv = L[[13]],
		sigma = L[[14]],
		semu = L[[15]],
		sesig = L[[16]],
		muth = L[[17]],
		varth = L[[18]],
		eff = L[[19]],
		doci = L[[20]],
		tol = L[[21]],
		guess = L[[22]],
		beta = L[[23]])
}
"setup"<-
function(lo, hi, adtim)
{
# This function is used within sis1 and gibbs1.
# Users need not worry about this function.
	infty <- -99
	cc <- c(lo, hi)
	last.inf <- any(cc == infty)
	cc <- unique(sort(cc[cc != infty]))
	time <- unique(sort(c(cc, adtim)))
	ntim <- length(time)
	nadtim <- length(time) - length(cc)
	intim <- match(cc, time)
	time <- c(time, infty)
	if(last.inf) {
		cc <- c(cc, infty)
		intim <- c(intim, ntim + 1)
	}
	inlo <- match(lo, cc)
	inhi <- match(hi, cc) - 1
	list(inlo = inlo, inhi = inhi, cc = cc, time = time, ntim = ntim, intim
		 = intim, nadtim = nadtim)
}
"sis1"<-
function(...)
{
	sisx("sis1", ...)
}
"sis1.call"<-
function(...)
{
	sisx.call("sis1", ...)
}
"sis1m"<-
function(data, prior, nrep, proposal = prior, seed = NULL, eps = 0.01, adtim = 
	NULL, ord = T, domcse = T, ci = 0, npass.ci = 3, nrep.ci = NULL, tol = 
	0, reuse = F)
{
#
# Required arguments:
# ------------------
# data is a matrix containing the data intervals.  The matrix has two
#   columns.  Each row of the matrix gives one data interval.  The time
#   value infinity may be coded as NA, Inf, or -99.
#   If column 2 of data consists entirely of 1's and 0's  we assume
#    these are indicators of death and right-censoring respectively.
# prior is a vector specifying the prior distribution.  It contains
#   three values a, b, totm.  The values a, b specify the gamma mixing
#   distribution.  The value totm is the total mass alpha(R) of the 
#   Dirichlet measure.
# nrep is the number of iterations of the SIS sampler to record.
#
# Optional arguments:
# ------------------
# proposal is a vector of three values specifying the proposal 
#   distribution used in the sequential importance sampling.  The 
#   values are similar to those in the argument prior.
# seed is a vector of three integer seeds for the random number 
#    generators.  The first seed should be between 1 and 32.
#    If no seeds are given, some are chosen at random.
# eps is a small positive value.  The program currently does not
#    handle exact observed death times t, but replaces these by 
#    interval censored observations (t-eps,t).
# ord is a logical variable (T or F).  When ord = T, the data is
#    ordered according to the prior probability of the data intervals.
#    When ord = F, the data is not ordered.
# adtim is a vector of additional time points at which the survival
#    function is to be computed.  If this is not specified, then
#    the survival function is only computed at the endpoints of the
#    data intervals.
#
	infty <- -99
	tmp <- lohi(data, eps)
	lo <- tmp$lo
	hi <- tmp$hi
	nn <- length(lo)
	aa0 <- prior[1]
	bb0 <- prior[2]
	totm0 <- prior[3]
	aa1 <- proposal[1]
	bb1 <- proposal[2]
	totm1 <- proposal[3]
	if(ord) {
		phi <- bb0^aa0/(bb0 + hi)^aa0
		phi[hi == infty] <- 0
		plo <- bb0^aa0/(bb0 + lo)^aa0
		prb <- plo - phi
		prb[tmp$tie] <- 0
		pi <- order(prb)
		lo <- lo[pi]
		hi <- hi[pi]
	}
	start.ranf(seed = seed)
	misc <- setup(lo, hi, adtim)
	inlo <- misc$inlo
	inhi <- misc$inhi
	cc <- misc$cc
	nc <- length(cc)
	intim <- misc$intim
	time <- misc$time
	ntim <- misc$ntim
	surv <- rep(0, ntim)
	sigma <- rep(0, ntim)
	if(domcse)
		semu <- rep(0, ntim)
	else semu <- 0
	ess <- 0
	pmax <- 0
	muth <- 0
	varth <- 0
	L <- list(lo = lo, hi = hi, inlo = inlo, inhi = inhi, nn = nn, aa0 = 
		aa0, bb0 = bb0, totm0 = totm0, aa1 = aa1, bb1 = bb1, totm1 = 
		totm1, nrep = nrep, cc = cc, intim = intim, nc = nc, domcse = 
		domcse, time = time, ntim = ntim, surv = surv, sigma = sigma, 
		semu = semu, ess = ess, pmax = pmax, muth = muth, varth = varth,
		doci = F, tol = tol, guess = 0, beta = 0)
	ans <- sis1m.call(L)
	if(ci > 0)
		ans <- c(ans, ci = ci, cisnrxx(ans[1:25], ci, "sis1m.call", 
			nrep.ci, npass.ci, min.size = 500, tol = tol, reuse = 
			reuse))
	class(ans) <- "mixdir"
	ans
}
"sis1m.call"<-
function(L)
{
	storage.mode(L[[1]]) <- "double"
	storage.mode(L[[2]]) <- "double"
	storage.mode(L[[3]]) <- "integer"
	storage.mode(L[[4]]) <- "integer"
	storage.mode(L[[5]]) <- "integer"
	storage.mode(L[[6]]) <- "double"
	storage.mode(L[[7]]) <- "double"
	storage.mode(L[[8]]) <- "double"
	storage.mode(L[[9]]) <- "double"
	storage.mode(L[[10]]) <- "double"
	storage.mode(L[[11]]) <- "double"
	storage.mode(L[[12]]) <- "integer"
	storage.mode(L[[13]]) <- "double"
	storage.mode(L[[14]]) <- "integer"
	storage.mode(L[[15]]) <- "integer"
	storage.mode(L[[16]]) <- "logical"
	storage.mode(L[[17]]) <- "double"
	storage.mode(L[[18]]) <- "integer"
	storage.mode(L[[19]]) <- "double"
	storage.mode(L[[20]]) <- "double"
	storage.mode(L[[21]]) <- "double"
	storage.mode(L[[22]]) <- "double"
	storage.mode(L[[23]]) <- "double"
	storage.mode(L[[24]]) <- "double"
	storage.mode(L[[25]]) <- "double"
	storage.mode(L[[26]]) <- "logical"
	storage.mode(L[[27]]) <- "double"
	storage.mode(L[[28]]) <- "double"
	storage.mode(L[[29]]) <- "double"
	.Fortran("sis1m",
		lo = L[[1]],
		hi = L[[2]],
		inlo = L[[3]],
		inhi = L[[4]],
		nn = L[[5]],
		aa0 = L[[6]],
		bb0 = L[[7]],
		totm0 = L[[8]],
		aa1 = L[[9]],
		bb1 = L[[10]],
		totm1 = L[[11]],
		nrep = L[[12]],
		cc = L[[13]],
		intim = L[[14]],
		nc = L[[15]],
		domcse = L[[16]],
		time = L[[17]],
		ntim = L[[18]],
		surv = L[[19]],
		sigma = L[[20]],
		semu = L[[21]],
		ess = L[[22]],
		pmax = L[[23]],
		muth = L[[24]],
		varth = L[[25]],
		doci = L[[26]],
		tol = L[[27]],
		guess = L[[28]],
		beta = L[[29]])
}
"sis2"<-
function(...)
{
	sisx("sis2", ...)
}
"sis2.call"<-
function(...)
{
	sisx.call("sis2", ...)
}
"sisx"<-
function(algo, data, prior, nrep, seed = NULL, eps = 0.01, adtim = NULL, ord = 
	T, domcse = T, ci = 0, npass.ci = 3, nrep.ci = NULL, tol = 0, reuse = F
	)
{
#
# Required arguments:
# ------------------
# data is a matrix containing the data intervals.  The matrix has two
#   columns.  Each row of the matrix gives one data interval.  The time
#   value infinity may be coded as NA, Inf, or -99.
#   If column 2 of data consists entirely of 1's and 0's  we assume
#    these are indicators of death and right-censoring respectively.
# prior is a vector specifying the prior distribution.  It contains
#   three values a, b, totm.  The values a, b specify the gamma mixing
#   distribution.  The value totm is the total mass alpha(R) of the 
#   Dirichlet measure.
# nrep is the number of iterations of the SIS sampler to record.
#
# Optional arguments:
# ------------------
# seed is a vector of three integer seeds for the random number 
#    generators.  The first seed should be between 1 and 32.
#    If no seeds are given, some are chosen at random.
# eps is a small positive value.  The program currently does not
#    handle exact observed death times t, but replaces these by 
#    interval censored observations (t-eps,t).
# ord is a logical variable (T or F).  When ord = T, the data is
#    ordered according to the prior probability of the data intervals.
#    When ord = F, the data is not ordered.
# adtim is a vector of additional time points at which the survival
#    function is to be computed.  If this is not specified, then
#    the survival function is only computed at the endpoints of the
#    data intervals.
#
	infty <- -99
	tmp <- lohi(data, eps)
	lo <- tmp$lo
	hi <- tmp$hi
	nn <- length(lo)
	aa0 <- prior[1]
	bb0 <- prior[2]
	totm0 <- prior[3]
	if(ord) {
		phi <- bb0^aa0/(bb0 + hi)^aa0
		phi[hi == infty] <- 0
		plo <- bb0^aa0/(bb0 + lo)^aa0
		prb <- plo - phi
		prb[tmp$tie] <- 0
		pi <- order(prb)
		lo <- lo[pi]
		hi <- hi[pi]
	}
	start.ranf(seed = seed)
	misc <- setup(lo, hi, adtim)
	inlo <- misc$inlo
	inhi <- misc$inhi
	cc <- misc$cc
	nc <- length(cc)
	intim <- misc$intim
	time <- misc$time
	ntim <- misc$ntim
	surv <- rep(0, ntim)
	sigma <- rep(0, ntim)
	if(domcse)
		semu <- rep(0, ntim)
	else semu <- 0
	ess <- 0
	pmax <- 0
	muth <- 0
	varth <- 0
	L <- list(lo = lo, hi = hi, inlo = inlo, inhi = inhi, nn = nn, aa0 = 
		aa0, bb0 = bb0, totm0 = totm0, nrep = nrep, cc = cc, intim = 
		intim, nc = nc, domcse = domcse, time = time, ntim = ntim, surv
		 = surv, sigma = sigma, semu = semu, ess = ess, pmax = pmax, 
		muth = muth, varth = varth, doci = F, tol = tol, guess = 0, 
		beta = 0)
	func <- paste(algo, "call", sep = ".")
	ans <- eval(call(func, L))
	if(ci > 0)
		ans <- c(ans, ci = ci, cisnrxx(ans[1:22], ci, func, nrep.ci,
			npass.ci, min.size = 500, tol = tol, reuse = reuse))
	class(ans) <- "mixdir"
	ans
}
"sisx.call"<-
function(algo, L)
{
	storage.mode(L[[1]]) <- "double"
	storage.mode(L[[2]]) <- "double"
	storage.mode(L[[3]]) <- "integer"
	storage.mode(L[[4]]) <- "integer"
	storage.mode(L[[5]]) <- "integer"
	storage.mode(L[[6]]) <- "double"
	storage.mode(L[[7]]) <- "double"
	storage.mode(L[[8]]) <- "double"
	storage.mode(L[[9]]) <- "integer"
	storage.mode(L[[10]]) <- "double"
	storage.mode(L[[11]]) <- "integer"
	storage.mode(L[[12]]) <- "integer"
	storage.mode(L[[13]]) <- "logical"
	storage.mode(L[[14]]) <- "double"
	storage.mode(L[[15]]) <- "integer"
	storage.mode(L[[16]]) <- "double"
	storage.mode(L[[17]]) <- "double"
	storage.mode(L[[18]]) <- "double"
	storage.mode(L[[19]]) <- "double"
	storage.mode(L[[20]]) <- "double"
	storage.mode(L[[21]]) <- "double"
	storage.mode(L[[22]]) <- "double"
	storage.mode(L[[23]]) <- "logical"
	storage.mode(L[[24]]) <- "double"
	storage.mode(L[[25]]) <- "double"
	storage.mode(L[[26]]) <- "double"
	.Fortran(algo,
		lo = L[[1]],
		hi = L[[2]],
		inlo = L[[3]],
		inhi = L[[4]],
		nn = L[[5]],
		aa0 = L[[6]],
		bb0 = L[[7]],
		totm0 = L[[8]],
		nrep = L[[9]],
		cc = L[[10]],
		intim = L[[11]],
		nc = L[[12]],
		domcse = L[[13]],
		time = L[[14]],
		ntim = L[[15]],
		surv = L[[16]],
		sigma = L[[17]],
		semu = L[[18]],
		ess = L[[19]],
		pmax = L[[20]],
		muth = L[[21]],
		varth = L[[22]],
		doci = L[[23]],
		tol = L[[24]],
		guess = L[[25]],
		beta = L[[26]])
}
"splice.matrices"<-
function(...)
{
	x <- list(...)
	ans <- c()
	for(i in x)
		ans <- cbind(ans, i)
	matrix(as.vector(t(ans)), ncol(x[[1]]), length(x) * nrow(x[[1]]))
}
"start.ranf"<-
function(seed = NULL)
{
# seed is a vector of 3 positive integers.
# seed[1] must be less than or equal to 32.
#------------------------------------------------------------
# First check to see if ranf has already been initialized.
	if(is.null(seed) && .Fortran("rgnqsd",
		T)[[1]]) return(NULL)
	if(is.null(seed)) {
		i1 <- ceiling(32 * runif(1))
		i2 <- ceiling(9999 * runif(1))
		i3 <- ceiling(9999 * runif(1))
	}
	else {
		i1 <- seed[1]
		if(i1 > 32)
			i1 <- (i1 %% 32) + 1
		i2 <- seed[2]
		i3 <- seed[3]
	}
	.Fortran("setall",
		as.integer(i2),
		as.integer(i3))
	.Fortran("setcgn",
		as.integer(i1))
	NULL
}
"summary.mixdir"<-
function(L, prn = T)
{
# Set prn=F to turn off automatic printing of ess and pmax.
	if(prn && (length(L$ess) > 0)) {
		print(paste("ess =", format(L$ess), " pmax =", format(L$pmax)), 
			quote = F)
		cat("\n")
	}
	ntim <- L$ntim
	ans <- cbind(time = L$time[1:ntim], mean = L$surv)
	if(length(L$semu) == ntim)
		ans <- cbind(ans, se.mean = L$semu)
	if(!is.null(lim <- L$limits))
		ans <- cbind(ans, lower = lim[1,  ], upper = lim[2,  ])
	if(length(L$sigma) == ntim)
		ans <- cbind(ans, sigma = L$sigma)
	if(length(L$sesig) == ntim)
		ans <- cbind(ans, se.sigma = L$sesig)
	ans
}
"ranf"<-function(n){
	.Fortran("vranf",as.integer(n),double(n))[[2]]
}