function( iters, filename)
{

#########################################
# data should be in list saved to file with dput
# x1 needs to include a column of 1's for intercept
# x2 must NOT include the column of 1's
datalist <- dget(filename)
y <- datalist$y
datalist <- dget(filename)
y1 <- datalist$s

n <- length(y1)
x1 <- cbind(rep(1,n), datalist$trt, datalist$strat)
x2 <- x1[ ,2:3]
status <- as.numeric( !(is.na(datalist$surv)))  # 1 = event,0 = censored
ncens <- n - sum(status)
y2 <- datalist$surv
y2[ status < 1 ] <- datalist$cent[status<1]
y <- cbind(y1,y2)
print(status)
print(y)

#x1 <- datalist$x1
#x2 <- datalist$x2
#status <- datalist$status	# 1 = event,0 = censored
#n <- nrow(y) 
#ncens <- n - sum(status)
#x2 <- matrix( x2, nrow = n)
print(c("nrow(y)", nrow(y), "ncol(y)", ncol(y)))
iters <- iters + 1		# allow initial values to be iter 1

library(survival)
survout<-survreg(Surv(y[,2], status) ~ x2[,1] + x2[,2])
# print(survout$linear.predictor)
print(summary(survout))
print(summary(survreg(Surv(y[,2], status) ~ x2[,1] + x2[,2], scale=1)))
print(summary(lm(y[,1] ~ x2[,1] + x2[,2])))




###########################################
R <- matrix( c(1,0,0,.1), nrow = 2)
nu <- 4
# parms of inverse Wishart prior on cov mat
# take flat prior on beta

R221 <- R[2,2] - R[1,2]^2 / R[1,1]
a1 <- (n + nu - 1) / 2
###########################################
# one-time calculations

p1 <- ncol(x1)  # number of predictors (including intercept) for y1
p2 <- ncol(x2)  # number of predictors (NOT including intercept) for y2

print( c("p1", p1, "p2", p2)) 
x1tx1inv <- solve( t(x1[,2:p1]) %*% x1[,2:p1] )
cholx1tx1inv <- chol(x1tx1inv)

y2copy <- y[,2]	# keep values of censoring times in y2copy ; impute into y[,2]

#x2tx2inv <- solve( t(x2[,2:p2]) %*% x2[,2:p2] )
#cholx2tx2inv <- chol(x2tx2inv)
###########################################

# define functions used in computing acceptance probabilities

###########
# loglik1 #
###########
loglik1 <- function( mu1i, mu2i, delta1, delta3, multi, y, alph)
{
# terms including beta1
 - (sum(y[,1] - mu1i)^2) / (2 * delta1) - sum(((y[,2]^alph - multi * mu2i)^2) / (2 * delta3 * multi ^ 2)) - sum( log(multi) ) - sum( log(1 - pnorm( -mu2i/sqrt(delta3) ) ))

# print(c("sum",sum( log(1 - pnorm( -mu2i/sqrt(delta3) ) ))))

#- sum( log(1 - pmin(pnorm( -mu2i/sqrt(delta3) ), rep(0.99999999, nrow(y)) )))

}
########### line 50
# loglik2 #
###########
loglik2 <- function( mu1i, mu2i,  delta3,  multi, y, alph)
{
# terms including beta2 and delta2 
- sum(((y[,2] ^ alph - multi * mu2i)^2) / (2 * delta3 * multi^2)) - sum(log(multi)) - sum( log(1 - pnorm( -mu2i/sqrt(delta3) ) ))
}

###########
# loglik3 #
###########
loglik3 <- function( mu1i, mu2i, delta3,  multi, y, alph)
{
# terms including delta3
-(n/2) * log(delta3) - sum(((y[,2]^ alph - multi * mu2i)^2) / (2 * multi^2 * delta3))  - sum( log(1 - pnorm( -mu2i/sqrt(delta3) ) ))
}

##############
# loglikalph #
##############
loglikalph <- function( mu2i, delta3,  multi, y, alph)
{
# terms including alph
nrow(y) * log(alph) + (alph-1) * sum(log(y[,2])) - sum(((y[,2]^ alph - multi * mu2i)^2) / (2 * multi^2 * delta3))
}

############
# loglikli #
############
loglikli <- function( i, mu1i, mu2i,  delta3,  multii, y, alph )
{
# terms involving lambda_i
- ((y[i,2]^alph - multii * mu2i[i])^2) / (2 * delta3 * multii^2) - log(multii) 
}

#############
# logprior1 #
#############
logprior1 <- function( R, nu, d1)
{
-( (nu - 1) /  2 + 1) * log(d1) - R[1,1] / (2 * d1)
}
#############
# logprior2 #
#############
logprior2 <- function( R, nu, d2, d3)
{
- log(d3) / 2 - (R[1,1] / (2 * d3) ) * (d2 - (R[1,2]/R[1,1]))^2
}
############# line 100
# logprior3 #
#############
logprior3 <- function( R, R221, nu, d2, d3)
{
-( (nu + 1) /  2 + 1) * log(d3)  - (R[1,1]  * (d2 - (R[1,2]/R[1,1]))^2 
+ R221) / (2 * d3) 
}

#############
# logpriorli #
#############
logpriorli <- function(li)
{
# exponential(1)

-li
}

###############
# logprioralph #
###############
logprioralph <- function(alph)
{
# exponential(0.001)

- 0.001 * alph
}


#################
# logtruncdenom #
#################
logtruncdenom <- function(mu1i, mu2i, delta3)
{
-sum( log(1 - pnorm( -mu2i/sqrt(delta3) ) )) # added minus sign 07/28
}



############################################
# initialize output matrices

beta1 <- matrix(rep(0, (p1) * iters), nrow = iters)
beta2 <- matrix(rep(0.5, (p2) * iters), nrow=iters)
delta <- matrix( c(rep(1,iters), rep(0,iters), rep(200,iters)), nrow=iters)
beta1[1,] <- c(0.177, 0.356, 0.191)
beta2[1,] <- c(-0.7, -0.95)

Sigma <- matrix( rep(0, 3 * iters), nrow = iters)
lambda <- matrix(rep(1, n * iters), nrow = iters)
beta20 <- matrix( rep(0,iters), ncol = 1)
correl <- matrix(rep(0,iters), ncol=1)
alph <- matrix(rep(1.34,iters), ncol = 1)
#############################################
# iterate
# line 150
for (i in 2:iters) {

	# generate exact values of censored y2's; must be greater than censoring times

	mu1i <- x1 %*% beta1[i-1,]
	mu2i <- delta[i-1,2] * (y[,1] - mu1i)
	multi <- sqrt(2 * lambda[i-1,]) * exp( - x2 %*% beta2[i-1, ])

	y[ status < 1, 2] <- (gentruncnorm0807( ncens, (mu2i[status < 1] * multi[status < 1]), 
multi[status < 1] * sqrt(delta[i-1,3]), y2copy[status < 1]^alph[i-1], 10000, 0, 1)) ^ ( 1 / alph[i-1])

# print(cbind(y[status<1,2], y2copy[status<1]))
	# generate delta1 from its inverse gamma full condl

	b1 <-  (( sum((y[,1] - mu1i)^2) + R[1,1]) / 2)
	delta[i,1] <- 1 / rgamma( 1, a1, b1)

	print(c("delta[i,1]", delta[i,1], "a1", a1, "b1", b1))
	##################################

# generate delta3   
	cand <- exp( rnorm(1, log(delta[i-1,3]), 0.1) )
	# log normal is not symmetric; so cand-gen density will
	# be part of acceptance prob
	u <- runif(1)

	logalpha <- loglik3( mu1i, mu2i, cand,  multi, y, alph[i-1])- loglik3( mu1i, mu2i, delta[i-1,3], multi, y, alph[i-1])

	logalpha <- logalpha + logprior3(R, R221, nu, delta[i-1,2], cand) - logprior3(R, R221, nu, delta[i-1,2], delta[i-1,3] )
	logalpha <- logalpha + log(cand) - log(delta[i-1,3])
	if (log(u) < logalpha )
		delta[i,3] <- cand
		
	else
		{delta[i,3] <- delta[i-1,3]
	
		}
print(c("delta[i,3]", delta[i,3]))
	##################################
# generate delta2
	# cand gen density: normal centered at previous value

	cand <- rnorm( 1, delta[i-1,2], 0.1)
	mu2inew <- mu2i + (y[,1] - mu1i) * (cand - delta[i-1,2])
	logalpha <- loglik2( mu1i, mu2inew, delta[i,3],  multi, y, alph[i-1]) - loglik2( mu1i, mu2i,  delta[i,3],  multi, y, alph[i-1])
	logalpha <- logalpha + logprior2( R, nu, cand, delta[i,3]) - logprior2( R, nu, delta[i-1,2], delta[i,3])
	u <- runif(1)
	if (log(u) < logalpha )
		{
		delta[i,2] <- cand
		mu2i <- mu2inew
		}
	else
		delta[i,2] <- delta[i-1,2]
			
print(c("delta[i,2]", delta[i,2]))
	##################################
# generate beta1
	# first the intercept
	# cand-gen density: normal centered at old value
	cand <- rnorm( 1, beta1[i-1,1], 0.1 )
	mu1inew <- mu1i - beta1[i-1,1] + cand
	mu2inew <- mu2i + delta[i,2] * (mu1i - mu1inew)
print("here")
	logalpha <- loglik1(mu1inew, mu2inew, delta[i,1], delta[i,3], multi, y, alph[i-1]) - loglik1( mu1i, mu2i, delta[i,1], delta[i,3], multi, y, alph[i-1] )
print(loglik1(mu1inew, mu2inew, delta[i,1], delta[i,3], multi, y, alph[i-1]))
print(loglik1( mu1i, mu2i, delta[i,1], delta[i,3], multi, y, alph[i-1] ))
	u <- runif(1)
 	if (log(u) < logalpha){
		beta1[i,1] <- cand
		mu1i <- mu1inew
		mu2i <- mu2inew
		}
	else
		beta1[i,1] <- beta1[i-1,1]
print(c("beta10", beta1[i,1]))
	# now the slopes
	# cand gen density: full condl if there were no truncation
	# this is normal, mean betahat, cov mat delta1 * x1tx1inv

	betahat <- x1tx1inv %*% t(x1[,2:p1]) %*% (y[,1] - beta1[i,1])
	z <- rnorm( (p1-1), 0, sd=sqrt(delta[i,1]) )
	cand <- betahat + t(cholx1tx1inv) %*% z
	mu1inew <-  x1 %*% c(beta1[i,1], cand)
	mu2inew <- mu2i + delta[i,2] * (mu1i - mu1inew)
	logalpha <- loglik2( mu1inew, mu2inew, delta[i,3],  multi, y, alph[i-1]) - loglik2( mu1i, mu2i,  delta[i,3],  multi, y, alph[i-1])
	u <- runif(1)
 	if (log(u) < logalpha){
		beta1[i,2:p1] <- cand
		mu1i <- mu1inew
		mu2i <- mu2inew
		}
	else
		beta1[i,2:p1] <- beta1[i-1,2:p1]

print(c("betahat", betahat, "cand", cand, "logalpha", logalpha, beta1[i,2:p1]))
	##################################
# generate beta2
	# now the slopes
	# cand gen density: normal centered at old value

	cand <- rnorm( p2, beta2[i-1, ], 0.1)
	multinew <- sqrt(2 * lambda[i-1,]) * exp( - x2 %*% cand)	
	
	logalpha <- loglik2(mu1i, mu2i, delta[i,3], multinew, y, alph[i-1]) - loglik2(mu1i, mu2i, delta[i,3], multi, y, alph[i-1])
	u <- runif(1)
 	if (log(u) < logalpha){
		beta2[i,] <- cand
		multi <- multinew
		}
	else
		beta2[i,] <- beta2[i-1,]
print(c("beta2", beta2[i,]))
# generate lambda
	for (j in 1:n) {
		cand <- exp( rnorm(1, log(lambda[i-1,j]), 0.1) )
		multinewi <- sqrt(2 * cand) * exp( - x2[j,] %*% beta2[i,])
		logalpha <- loglikli( j, mu1i, mu2i,  delta[i,3],  multinewi, y, alph[i-1] ) - loglikli( j, mu1i, mu2i,  delta[i,3],  multi[j], y, alph[i-1] )
		logalpha <- logalpha + logpriorli(cand) - logpriorli(lambda[i-1,j])
		logalpha <- logalpha + log(cand) - log(lambda[i-1,j])
 		if (log(u) < logalpha){
			lambda[i,j] <- cand
			multi[j] <- multinewi
			}
		else
			lambda[i,j] <- lambda[i-1,j]
	}
print("done with lambdas")
# generate alph
	cand <- exp( rnorm(1, log(alph[i-1]), 0.1) )

		logalpha <- loglikalph(mu2i,  delta[i,3],  multi, y, cand ) - loglikalph(mu2i,  delta[i,3],  multi, y, alph[i-1] )
		logalpha <- logalpha + logprioralph(cand) - logprioralph(alph[i-1,])
		logalpha <- logalpha + log(cand) - log(alph[i-1])
 		if (log(u) < logalpha){
			alph[i] <- cand
			}
		else
			alph[i] <- alph[i-1]
print(c("alph", i, alph[i]))

Sigma[i,1] <- delta[i,1]
Sigma[i,2] <- delta[i,2] * delta[i,1]
Sigma[i,3] <- delta[i,3] + Sigma[i,2]^2 / Sigma[i,1]
beta20[i,1] <- log(Sigma[i,3]) / 2
correl[i,1] <- Sigma[i,2] / sqrt( Sigma[i,1] * Sigma[i,3])
}

#myobj <- list(y=round(y,4))
#dput(myobj)
	
output <- cbind(beta1, beta20, beta2, alph, delta, Sigma, correl)
print("beta2")
print(beta2)
print("beta20")
print(beta20)
print( apply( output[ ((iters/2) + 1):iters, ], 2, mean))
print( apply( output[ ((iters/2) + 1):iters, ], 2, sd))
#print(var( output[ ((iters/2) + 1):iters, ]))
output

}