# modified Bessel function of fractional order

bessel <- function(y, lambda)
{
	zz <- list()
	if(length(y)>1&&length(lambda)>1&&length(y)!=length(lambda))
		stop("Bessel: length of y must equal length of lambda")
	else if(length(y)==1&&length(lambda)>1)y <- rep(y,length(lambda))
	else if(length(y)>1&&length(lambda)==1)lambda <- rep(lambda,length(y))
	z <- .C("bessik",
		y=y,
		lambda=abs(lambda),
		n=as.integer(length(y)),
		ri=double(length(y)),
		rk=double(length(y)),
		rip=double(length(y)),
		rkp=double(length(y)),
		err=integer(1))
	if(z$err>0)warning("Bessel: non-zero error code")
	z$ri <- ifelse(lambda<0,z$ri+2*sin(abs(lambda)*pi)*z$rk/pi,z$ri)
	z
}
pinvgauss <- function(q, m, s){
	t <- q/m
	v <- sqrt(q*s)
	pnorm((t-1)/v)+exp(2/(m*s))*pnorm(-(t+1)/v)}

dinvgauss <- function(y, m, s)
	exp(-(y-m)^2/(y*s*m^2)/2)/sqrt(2*pi*s*y^3)

plaplace <- function(q, m=0, s=1){
	u <- (q-m)/s
	t <- exp(-abs(u))/2
	ifelse(u<0,t,1-t)}

dlaplace <- function(y, m=0, s=1) exp(-abs(y-m)/s)/(2*s)

ddoublepois <- function(y, m, s){
	if(any(y<0))stop("y must contain non-negative values")
	if(any(m<=0))stop("m must be positive")
	if(any(s<=0))stop("s must be positive")
	if(length(m)!=length(y)&&length(m)!=1)stop("m must be a scalar or a vector of the same length as y")
	if(length(s)!=length(y)&&length(s)!=1)stop("s must be a scalar or a vector of the same length as y")
	if(length(m)!=length(s))stop("m and s must have the same length")
	.C("ddp",
	as.integer(y),
	as.integer(3*max(c(y,100))),
	as.double(m),
	as.double(s),
	as.integer(length(y)),
	as.integer(length(s)>1),
	res=double(length(y)))$res}

dmultpois <- function(y, m, s){
	if(any(y<0))stop("y must contain non-negative values")
	if(any(m<=0))stop("m must be positive")
	if(any(s<=0))stop("s must be positive")
	if(length(m)!=length(y)&&length(m)!=1)stop("m must be a scalar or a vector of the same length as y")
	if(length(s)!=length(y)&&length(s)!=1)stop("s must be a scalar or a vector of the same length as y")
	if(length(m)!=length(s))stop("m and s must have the same length")
	.C("dmp",
	as.integer(y),
	as.integer(3*max(c(y,100))),
	as.double(m),
	as.double(s),
	as.integer(length(y)),
	as.integer(length(s)>1),
	res=double(length(y)))$res}

ddoublebinom <- function(y, m, s){
	if(is.vector(y)&&length(y)==2)y <- matrix(y,nrow=1)
	else if(!is.matrix(y)||ncol(y)!=2)stop("y must be a two column matrix")
	if(any(y<0))stop("y must contain non-negative values")
	if(length(m)!=nrow(y)&&length(m)!=1)stop("m must be a scalar or a vector of the same length as y")
	if(any(m<=0)||any(m>=1))stop("m must lie between 0 and 1")
	if(length(s)!=nrow(y)&&length(s)!=1)stop("s must be a scalar or a vector of the same length as y")
	if(any(s<=0))stop("s must be positive")
	if(length(m)!=length(s))stop("m and s must have the same length")
	.C("ddb",
	as.integer(y[,1]),
	as.integer(y[,2]),
	as.double(m),
	as.double(s),
	as.integer(nrow(y)),
	as.integer(length(s)>1),
	res=double(nrow(y)))$res}

dmultbinom <- function(y, m, s){
	if(is.vector(y)&&length(y)==2)y <- matrix(y,nrow=1)
	else if(!is.matrix(y)||ncol(y)!=2)stop("y must be a two column matrix")
	if(any(y<0))stop("y must contain non-negative values")
	if(length(m)!=nrow(y)&&length(m)!=1)stop("m must be a scalar or a vector of the same length as y")
	if(any(m<=0)||any(m>=1))stop("m must lie between 0 and 1")
	if(length(s)!=nrow(y)&&length(s)!=1)stop("s must be a scalar or a vector of the same length as y")
	if(any(s<=0))stop("s must be positive")
	if(length(m)!=length(s))stop("m and s must have the same length")
	.C("dmb",
	as.integer(y[,1]),
	as.integer(y[,2]),
	as.double(m),
	as.double(s),
	as.integer(nrow(y)),
	as.integer(length(s)>1),
	res=double(nrow(y)))$res}
fit.dist <- function(y, ni, dist="normal", breaks=F, delta=1, censor=F,
	exact=T, plot=F, add=F, main, xlab, ...)
{
	if(!missing(dist))dist <- match.arg(dist,c("binomial","Poisson",
		"negative binomial","geometric","zeta","normal","log normal",
		"inverse Gauss","logistic","exponential","Pareto","gamma",
		"Weibull"))
	if(missing(main))main <- paste("Histogram of", deparse(substitute(y)))
	if(missing(xlab))xlab <- paste(deparse(substitute(y)))
	n <- length(ni)
	if (length(delta)==1) delta <- rep(delta,n)
	if (breaks){
		if(length(y)!=n+1)
			stop("Breaks vector must be one longer than frequency vector")
		yi <- (y[1:n]+y[2:(n+1)])/2
		delta <- diff(y)}
	else yi <- y
	pi.hat <- ni/sum(ni)
        ybar <- weighted.mean(yi,ni)
	s2 <- weighted.mean((yi-ybar)^2,ni)
	if (dist=="binomial"){
		m <- length(yi)-1
		nu <- ybar/m
		pi.tilde <- gamma(m+1)/gamma(yi+1)/gamma(m-yi+1)*
			nu^yi*(1-nu)^(m-yi)
		param <- nu
		names(param) <- "nu.hat"
		p <- 1}
	else if (dist=="Poisson"){
		pi.tilde <- exp(-ybar)*ybar^yi/gamma(yi+1)
		param <- ybar
		names(param) <- "mu.hat"
		p <- 1}
	else if (dist=="geometric"){
		nu <- 1/(1+ybar)
		pi.tilde <- nu*(1-nu)^yi
		param <- nu
		names(param) <- "nu.hat"
		p <- 1}
	else if (dist=="negative binomial"){
		nu <- ybar/s2
		gam <- ybar^2/(s2-ybar)
		if(exact){
			nu <- log(nu/(1-nu))
			fcn <- function(p)
				-sum(ni*(lgamma(yi+p[2])-lgamma(p[2])+
				p[2]*p[1]-(yi+p[2])*log(1+exp(p[1]))))
#				p[2]*log(p[1])+yi*log(1-p[1])))
			z <- nlm(fcn, p=c(nu, gam), stepmax=sqrt(nu^2+gam^2)/2)
			nu <- exp(z$estimate[1])
			nu <- nu/(1+nu)
			gam <- z$estimate[2]}
		pi.tilde <- gamma(yi+gam)/gamma(yi+1)/gamma(gam)*
			nu^gam*(1-nu)^yi
		param <- c(nu, gam)
		names(param) <- c("nu.hat","gamma.hat")
		p <- 2}
	else if (dist=="zeta"){
		pi.tilde <- 1/yi
		nu <- sum(pi.tilde)
		pi.tilde <- pi.tilde/nu
		rho <- round(pi.hat[1]/pi.tilde[1]+0.1)
		if(exact){
			fcn <- function(p) {
				if(censor) const <- sum(yi^(-p[1]))
				else const <- sum(seq(1,30)^(-p[1]))
				sum(ni*(p[1]*log(yi)+log(const)))}
			z <- nlm(fcn, p=rho, stepmax=1)
			rho <- z$estimate[1]}
		pi.tilde <- yi^(-rho)
		nu <- sum(pi.tilde)
		pi.tilde <- pi.tilde/nu
		param <- rho
		names(param) <- "rho.hat"
		p <- 1}
        else if (dist=="normal"){
		mu.hat <- ybar
		sigma2.hat <- s2
		pi.tilde <- exp(-(yi-mu.hat)^2/(2*sigma2.hat))/
			sqrt(2*pi*sigma2.hat)
		param <- c(mu.hat,sigma2.hat)
		names(param) <- c("mu.hat","sigma2.hat")
		p <- 2}
        else if (dist=="log normal"){
		mu.hat <- weighted.mean(log(yi),ni)
		sigma2.hat <- weighted.mean((log(yi)-mu.hat)^2,ni)
		pi.tilde <- exp(-(log(yi)-mu.hat)^2/(2*sigma2.hat))/
			(yi*sqrt(2*pi*sigma2.hat))
		param <- c(mu.hat,sigma2.hat)
		names(param) <- c("mu.hat","sigma2.hat")
		p <- 2}
        else if (dist=="inverse Gauss"){
		mu.hat <- ybar
		sigma2.hat <- weighted.mean(1/yi,ni)-(1/ybar)
		pi.tilde <- exp(-(yi-mu.hat)^2/(2*yi*sigma2.hat*mu.hat^2))/
			sqrt(2*pi*yi^3*sigma2.hat)
		param <- c(mu.hat,sigma2.hat)
		names(param) <- c("mu.hat","sigma2.hat")
		p <- 2}
	else if (dist=="logistic"){
		mu.hat <- ybar
		sigma <- sqrt(s2)
		if(exact){
			fcn <- function(p)
				sum(ni*(pi*(yi-p[1])/p[2]/sqrt(3)+log(p[2])
				+2*log(1+exp(-pi*(yi-p[1])/p[2]/sqrt(3)))))
			z <- nlm(fcn, p=c(mu.hat,sigma), stepmax=10)
			mu.hat <- z$estimate[1]
			sigma <- z$estimate[2]}
		pi.tilde <- exp(-pi*(yi-mu.hat)/(sigma*sqrt(3)))
		pi.tilde <- pi*pi.tilde/(sigma*sqrt(3)*(1+pi.tilde)^2)
		param <- c(mu.hat,sigma)
		names(param) <- c("mu.hat","sigma.hat")
		p <- 2}
        else if (dist=="exponential"){
		pi.tilde <- (1/ybar)*exp(-yi/ybar)
		param <- ybar
		names(param) <- "mu.hat"
		p <- 1}
        else if (dist=="Pareto"){
		delta.hat <- yi[1]-delta[1]/2
		alpha.hat <- sum(ni)/sum(ni*log(yi/delta.hat))
		pi.tilde <- alpha.hat*delta.hat^alpha.hat/(yi^(alpha.hat+1))
		param <- c(alpha.hat,delta.hat)
		names(param) <- c("alpha.hat","delta.hat")
		p <- 2}
        else if (dist=="gamma"){
		alpha.hat <- ybar^2/s2
		mu.hat <- ybar
		if(exact){
			fcn <- function(p)
				-sum(ni*(p[2]*log(p[2])
				+(p[2]-1)*log(yi)-p[2]*log(p[1])
				-p[2]*yi/p[1]-lgamma(p[2])))
			z <- nlm(fcn, p=c(mu.hat,alpha.hat), stepmax=10)
			mu.hat <- z$estimate[1]
			alpha.hat <- z$estimate[2]}
		pi.tilde <- (alpha.hat^alpha.hat)*(yi^(alpha.hat-1))/
			(mu.hat^alpha.hat)*exp(-alpha.hat*yi/mu.hat)/
			gamma(alpha.hat)
		param <- c(alpha.hat,mu.hat)
		names(param) <- c("alpha.hat","mu.hat")
		p <- 2}
        else if (dist=="Weibull"){
		tamp <- ybar^2/(s2+ybar^2)
		Alpha.Weibull.fn <- function(y){
			alpha.trans.fn <- function(al)
				gamma(1+1/al)*gamma(1+1/al)/gamma(1+2/al)
			tol <- 0.001
			al.start <- 0.0001
			al.end <- 50
			al.mid <- 0.5*(al.start+al.end)
			y.tamp <- alpha.trans.fn(al.mid)
			while (abs(y.tamp-y)>tol){
				if ((y.tamp-y)>0) al.end <- al.mid
				else al.start <- al.mid
				al.mid <- 0.5*(al.start+al.end)
				y.tamp <- alpha.trans.fn(al.mid)}
			al.mid}
		alpha.hat <- Alpha.Weibull.fn(tamp)
		mu.hat <- weighted.mean(yi^alpha.hat,ni)
		if(exact){
			fcn <- function(p)
				-sum(ni*(log(p[2])+(p[2]-1)*log(yi)
				-log(p[1])-yi^p[2]/p[1]))
			z <- nlm(fcn, p=c(mu.hat,alpha.hat), stepmax=10)
			mu.hat <- z$estimate[1]
			alpha.hat <- z$estimate[2]}
		pi.tilde <- alpha.hat*yi^(alpha.hat-1)/mu.hat*
			exp(-yi^alpha.hat/mu.hat)
		param <- c(alpha.hat,mu.hat)
		names(param) <- c("alpha.hat","mu.hat")
		p <- 2}
	pi.tilde <- pi.tilde*delta
        if (censor)
		pi.tilde[length(pi.tilde)] <- 1-sum(pi.tilde[1:(length(pi.tilde)-1)])
	dev.comp <- rep(0,length(pi.tilde))
	dev.comp[ni>0] <- -2*ni[ni>0]*log(pi.tilde[ni>0]/pi.hat[ni>0])
	deviance <- sum(dev.comp)
	AIC <- deviance+2*p
	resid <- (ni-sum(ni)*pi.tilde)/sqrt(sum(ni)*pi.tilde)
	result.output <- c(ybar,s2,param,deviance,AIC)
	names(result.output) <- c("mean","variance",names(param),"deviance","AIC")
	cat(dist," distribution,","  n = ",sum(ni),"\n",sep="")
	print(result.output)
        cat("\n")
	if(censor)nn <- n-1
	else nn <- n
	if (plot){
		if (!add){
			cum.histo (c(yi-delta/2, yi[n]+delta[n]/2), ni, prob=T,
				ylab="Probability", main, xlab, ...)
			lines(yi[1:nn], pi.tilde[1:nn]/delta[1:nn])}
		else lines(yi[1:nn], pi.tilde[1:nn]/delta[1:nn], add=T, ...)}
	if(length(unique(delta))==1)
		df <- data.frame(yi,ni,pi.hat,pi.tilde,dev.comp,resid)
	else
		df <- data.frame(yi,ni,delta,pi.hat,pi.tilde,dev.comp,resid)
#	list(parameters=result.output,df=df)
	df}

cum.histo <- function (breaks, freq, prob = F,
	main = paste("Histogram of", deparse(substitute(breaks))),
	xlab = deparse(substitute(breaks)), ylab,
        xlim = range(breaks), ...)
{
	if (prob) {
                freq <- freq/(sum(freq) * diff(breaks))
                if (missing(ylab)) 
                        ylab <- "Relative Frequency"
        }
        else if (missing(ylab)) 
                ylab <- "Frequency"
	plot(breaks, c(freq, 0), type = "n", main = main,
		xlab = xlab, ylab = ylab, box = F, ...)
	rect(breaks[-length(breaks)], 0, breaks[-1], freq, border = par("fg"))}
fmr <- function(y, dist="normal", pmu=NULL, pmix=NULL, pshape=NULL, mu=NULL,
	mix=NULL, linear=NULL, censor="right",exact=F, wt=1, delta=1,
	print.level=0, typsiz=abs(p), ndigit=10, gradtol=0.00001,
	stepmax=10*sqrt(p%*%p), steptol=0.00001, iterlim=100, fscale=1){
dpnc <- function(y,m,s){
	r <- rep(0,length(s))
	for(i in 1:y)r <- r+exp(i*(1-s)*log(i)+i*s*log(m)+i*(s-1)-s*m-
		lgamma(i+1))
	r <- r+exp(-s*m)
	ifelse(wt,r,1)}
mpnc <- function(y,m,t){
	r <- rep(0,length(t))
	for(i in 0:y)r <- r+exp(i*log(m)+i*i*t-m-lgamma(i+1))
	ifelse(wt,r,1)}
fdbnc <- function(n,m,s){
	y <- 0:n
	a <- lchoose(n,y)+((n-y)*s)*log(1-m)+y*s*log(m)-n*(1-s)*log(n)
	y <- 1:n
	a <- a+c(0,y*(1-s)*log(y))
	y <- 0:(n-1)
	sum(exp(a+c((n-y)*(1-s)*log(n-y),0)))}
dbnc <- function(n,m,s){
	r <- rep(0,length(n))
	for(i in 1:length(n))r[i] <- fdbnc(n[i],m[i],s)
	ifelse(wt,r,1)}
mbnc <- function(n,m,s){
	r <- rep(0,length(n))
	for(i in 1:length(n))r[i] <- sum(exp(lchoose(n[i],0:n[i])+
		(n[i]-(0:n[i]))*log(1-m[i])+(0:n[i])*(log(m[i])+
		(n[i]-(0:n[i]))*s)))
	ifelse(wt,r,1)}
call <- sys.call()
if(!missing(dist)&&!is.function(dist)){
	dist <- match.arg(dist,c("binomial","beta binomial","double binomial",
	"mult binomial","Poisson","negative binomial","negative binomial",
	"mult Poisson","geometric","normal","inverse Gauss","logistic",
	"exponential","gamma","Weibull","extreme value","Cauchy",
	"Student t","Laplace"))}
if(!missing(pmu))npl <- length(pmu)
else npl <- 0
if(!missing(pmix))npm <- length(pmix)
else npm <- 0
if(dist!="binomial"&&dist!="Poisson"&&dist!="exponential"&&missing(pshape))
	stop("An estimate of the shape parameter must be given")
np <- npl+npm+1
p <- c(pmu,pmix,pshape)
if(is.function(dist)){
	fcn <- dist
	dist <- "own"}
if(is.list(y)&&!is.null(class(y))&&class(y)=="response"){
	if(is.null(y$censor))y <- y$y
	else y <- cbind(y$y,y$censor)}
if(dist=="Poisson"||dist=="negative binomial"||dist=="double Poisson"||dist=="mult Poisson"){
	if(!is.vector(y))stop("y must be a vector")
	n <- length(y)
	cens <- ifelse(y==0,1,0)}
else {
	if(length(dim(y))!=2||ncol(y)!=2)
		stop(paste("Two column matrix required for response:",
		if(dist=="binomial"||dist=="beta binomial"||dist=="double binomial"||dist=="mult binomial")"successes and failures"
		else "times and censor indicator"))
	else {
		n <- nrow(y)
		if(dist=="binomial"||dist=="beta binomial"||dist=="double binomial"||dist=="mult binomial"){
			if(missing(censor))
				stop("Censoring must be left, right, or both")
			if(censor!="left"&&censor!="right"&&censor!="both")
				stop("Censoring must be left, right, or both")
			lcens <- ifelse((censor=="left"|censor=="both")&
				y[,1]==0,1,0)
			rcens <-ifelse((censor=="right"|censor=="both")&
				y[,2]==0,1,0)
			if(censor=="both"){
				lcens <- lcens/2
				rcens <- rcens/2}
			nn <- y[,1]+y[,2]}
		else {
			if(any(delta<=0&y[,2]==1))
				stop("All deltas for uncensored data must be positive")
			else {
				delta <- ifelse(delta<=0,0.000001,delta)
				delta <- ifelse(y[,1]-delta/2<=0,delta-0.00001
				,delta)}
			y[,2] <- as.integer(y[,2])
			if(any(y[,2]!=-1&y[,2]!=0&y[,2]!=1))
				stop("Censor indicator must be -1, 0, or 1")
			if(censor!="left"&&censor!="right")
				stop("Censoring must be left or right")
			if(censor=="left"&!any(y[,2]==-1))
				stop("No left censored observations")
			if(censor=="right"&!any(y[,2]==0))
				stop("No right censored observations")
			cens <- as.integer(y[,2]==1)
			b <- as.integer((censor=="right"&y[,2]==0)|
				(censor=="left"&y[,2]==-1))
			r <- as.integer(censor=="left"&y[,2]==0)
			l <- as.integer(censor=="right"&y[,2]==-1)
			lc <- ifelse(censor=="left",1,0)
			rc <- ifelse(censor=="right",-1,1)}}
	if(dist=="double Poisson"||dist=="mult Poisson")
				my <- min(3*max(y),100)}
if((dist!="normal"&&dist!="logistic"&&dist!="Cauchy"&&
	dist!="Laplace"&&dist!="Student t"&&dist!="Poisson"&&
	dist!="negative binomial"&&dist!="double Poisson"&&
	dist!="mult Poisson")&&(any(y[,1]<=0)))
	stop("All response values must be > 0")
else if((dist=="Poisson"||dist=="negative binomial"||
		dist=="double Poisson"||dist=="mult Poisson")&&(any(y<0)))
		stop("All response values must be >= 0")
if(min(wt)<0)stop("All weights must be non-negative")
if(length(wt)==1)wt <- rep(wt,n)
if(length(delta)==1)delta <- rep(delta,n)
lin1 <- lin2 <- NULL
if(is.list(linear)){
	lin1 <- linear[[1]]
	lin2 <- linear[[2]]
	lmix <- is.language(lin2)}
else {
	lin1 <- linear
	lmix <- FALSE}
if(is.language(mu))lin1 <- mu
if(is.language(mix)){
	lin2 <- mix
	lmix <- is.language(lin2)}
nlp <- npl
if(is.language(lin1)){
	mt <- terms(lin1)
	if(is.numeric(mt[[2]])){
		dm1 <- matrix(1)
		colnames(dm1) <- "(Intercept)"
		npt1 <- 1
		if(!is.function(mu)){
			mu1 <- function(p) p[1]*rep(1,n)
			nlp <- 1}
		else mu1 <- function (p) mu(p, p[1]*rep(1,n))}
	else {
		mf <- model.frame(mt,sys.frame(sys.parent()),na.action=na.fail)
		dm1 <- model.matrix(mt, mf)
		npt1 <- ncol(dm1)
		if(!is.function(mu)){
			mu1 <- function(p) dm1 %*% p[1:npt1]
			nlp <- npt1}
		else mu1 <- function (p) mu(p, dm1 %*% p[1:npt1])}
	if(npl<npt1)stop("Not enough initial estimates for mu")}
else if(!is.function(mu)){
	mu1 <- function(p) p[1]*rep(1,n)
	nlp <- 1}
else mu1 <- mu
if(nlp!=npl)
	stop("Number of initial estimates for mu does not correspond to model")
npl1 <- npl+1
nlp <- npm
if(lmix){
	mt <- terms(lin2)
	if(is.numeric(mt[[2]])){
		dm2 <- matrix(1)
		colnames(dm2) <- "(Intercept)"
		npt2 <- 1
		if(!is.function(mix)) mixt <- function(p) {
			mf <- p[npl1]*rep(1,n)
			exp(mf)/(1+exp(mf))}
		else mixt <- function(p) {
			mf <- mix(p[npl1:np], p[npl1]*rep(1,n))
			exp(mf)/(1+exp(mf))}}
	else {
		mf <- model.frame(mt,sys.frame(sys.parent()),na.action=na.fail)
		dm2 <- model.matrix(mt, mf)
		npt2 <- ncol(dm2)
		if(!is.function(mix)){
			mixt <- function(p) {
				mf <- dm2 %*% p[npl1:(npl1+npt2-1)]
				exp(mf)/(1+exp(mf))}
			nlp <- npt2}
		else mixt <- function (p) {
			mf <- mix(p[npl1:np], dm2 %*% p[npl1:(nlp1+npt2-1)])
			exp(mf)/(1+exp(mf))}}
	if(npm<npt2)stop("Not enough initial estimates for mix")}
#else if(!is.function(mix)&&dist!="binomial"&&dist!="Poisson"&&
#	dist!="exponential"&&dist!="geometric"){
else if(!is.function(mix)){
	mixt <- function(p) exp(p[npl1])/(1+exp(p[npl1]))*rep(1,n)
	nlp <- 1}
else mixt <- function(p) {
	mf <- mix(p[npl1:np])
	exp(mf)/(1+exp(mf))}
if(nlp!=npm)
	stop("Number of initial estimates for mix does not correspond to model")
if(!is.numeric(mu1(pmu)))
	stop("The location function must return numerical values")
if(dist!="binomial"&&dist!="Poisson"&&dist!="exponential"&&
	dist!="geometric"&&!is.numeric(mixt(p)))
	stop("The mix function must return numerical values")
ret <- switch(dist,
	binomial={
		fcn <- function(p) {
			m <- mu1(p)
			s <- mixt(p)
			-wt*log((1-s)*(lcens+rcens)+s*m^y[,1]*(1-m)^y[,2])}
		const <- -wt*(lchoose(nn,y[,1]))},
	"beta binomial"={
		fcn <- function(p) {
			m <- mu1(p)
			s <- mixt(p)
			v <- exp(p[np])
			t <- v*m
			u <- v*(1-m)
			-wt*log((1-s)*(lcens+rcens)+s*beta(y[,1]+t,y[,2]+u)/
				beta(t,u))}
		const <- -wt*(lchoose(nn,y[,1]))},
	"double binomial"={
		fcn <- function(p) {
			-wt*log(.C("ddb",as.integer(y[,1]),as.integer(y[,2]),
				as.double(mu1(p)),as.double(exp(sh1(p))),
				as.integer(nrow(y)),as.integer(1),
				res=double(nrow(y)))$res)}
		const <- 0},
	"mult binomial"={
		fcn <- function(p) {
			-wt*log(.C("dmb",as.integer(y[,1]),as.integer(y[,2]),
				as.double(mu1(p)),as.double(exp(sh1(p))),
				as.integer(nrow(y)),as.integer(1),
				res=double(nrow(y)))$res)}
		const <- 0},
	Poisson={
		fcn <- function(p) {
			m <- mu1(p)
			s <- mixt(p)
			-wt*log((1-s)*cens+s*exp(-m)*m^y)}
		const <- wt*lgamma(y+1)},
	"negative binomial"={
		fcn <- function(p) {
			m <- mu1(p)
			s <- mixt(p)
			t <- exp(p[np])
			-wt*log((1-s)*cens+s*gamma(y+t)/gamma(t)
				*t^t*m^y/(t+m)^(y+t))}
		const <- wt*lgamma(y+1)},
	"double Poisson"={
		fcn <- function(p) {
			-wt*log(.C("ddp",as.integer(y),as.integer(my),
				as.double(mu1(p)),as.double(exp(sh1(p))),
				as.integer(length(y)),as.integer(1),
				res=double(length(y)))$res)}
		const <- 0},
	"mult Poisson"={
		fcn <- function(p) {
			-wt*log(.C("dmp",as.integer(y),as.integer(my),
				as.double(mu1(p)),as.double(exp(sh1(p))),
				as.integer(length(y)),as.integer(1),
				res=double(length(y)))$res)}
		const <- 0},
	normal={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- mixt(p)
				t <- exp(p[np]/2)
				pn <- pnorm(y[,1]-delta/2,m,t)
				-wt*log(s*cens*(pnorm(y[,1]+delta/2,m,t)-pn)
					+(1-cens)*((1+s*(rc*pn-lc))*b
					+s*(r+pn*(l-r))))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				s <- mixt(p)
				pn <- pnorm(y[,1]-delta/2,m,exp(p[np]/2))
				-wt*log(s*cens*exp((-(p[np]+(y[,1]-m)^2
					*exp(-p[np]))/2))
					+(1-cens)*((1+s*(rc*pn-lc))*b
					+s*(r+pn*(l-r))))}
			const <- wt*cens*(log(2*pi)/2-log(delta))}},
        "inverse Gauss"={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- mixt(p)
				t <- exp(p[np])
				pit <- pinvgauss(y[,1]-delta/2,m,t)
				-wt*log(s*cens*(pinvgauss(y[,1]+delta/2,m,t)-pit)
					+(1-cens)*((1+s*(rc*pit-lc))*b
					+s*(r+pit*(l-r))))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				s <- mixt(p)
				t <- exp(p[np])
				pit <- pinvgauss(y[,1]-delta/2,m,t)
				-wt*log(s*cens*exp(-(p[np]+(y[,1]-m)^2/
					(y[,1]*t*m^2))/2)
					+(1-cens)*((1+s*(rc*pit-lc))*b
					+s*(r+pit*(l-r))))}
			const <- wt*cens*(log(2*pi*y[,1]^3)/2-log(delta))}},
	logistic={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- mixt(p)
				t <- exp(p[np])*sqrt(3)/pi
				pl <- plogis(y[,1]-delta/2,m,t)
				-wt*log(s*cens*(plogis(y[,1]+delta/2,m,t)-pl)
					+(1-cens)*((1+s*(rc*pl-lc))*b
					+s*(r+pl*(l-r))))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				s <- mixt(p)
				t <- exp(p[np])*sqrt(3)/pi
				y1 <- (y[,1]-m)/t
				pl <- plogis(y[,1]-delta/2,m,t)
				-wt*log(s*cens*exp(-y1-log(t)
					-2*log(1+exp(-y1)))
					+(1-cens)*((1+s*(rc*pl-lc))*b
					+s*(r+pl*(l-r))))}
			const <- -wt*cens*log(delta)}},
        Student={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- mixt(p)
				t <- exp(p[np])
				ps <- pt(y[,1]-delta/2-m,t)
				-wt*log(s*cens*(pt(y[,1]+delta/2-m,t)-ps)
					+(1-cens)*((1+s*(rc*ps-lc))*b
					+s*(r+ps*(l-r))))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				s <- mixt(p)
				t <- exp(p[np])
				ps <- pt(y[,1]-delta/2-m,t)
				-wt*log(s*cens*gamma((t+1)/2)/gamma(t/2)*
					exp(-p[np]/2-((t+1)/2)*
					log(1+(y[,1]-m)^2/t))
					+(1-cens)*((1+s*(rc*ps-lc))*b
					+s*(r+ps*(l-r))))}
			const <- wt*cens*(log(pi)/2-log(delta))}},
	Cauchy={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- mixt(p)
				t <- exp(p[np]/2)
				pl <- pcauchy(y[,1]-delta/2,m,t)
				-wt*log(s*cens*(pcauchy(y[,1]+delta/2,m,t)-pl)
					+(1-cens)*((1+s*(rc*pl-lc))*b
					+s*(r+pl*(l-r))))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				s <- mixt(p)
				t <- exp(p[np]/2)
				pl <- pcauchy(y[,1]-delta/2,m,t)
				-wt*log(s*cens/(t*(1+(y[,1]-m)^2/t^2))
					+(1-cens)*((1+s*(rc*pl-lc))*b
					+s*(r+pl*(l-r))))}
			const <- -wt*cens*log(delta/pi)}},
        Laplace={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- mixt(p)
				t <- exp(p[np])
				pl <- plaplace((y[,1]-delta/2-m)/t)
				-wt*log(s*cens*(plaplace((y[,1]+delta/2-m)/t)-pl)
					+(1-cens)*((1+s*(rc*pl-lc))*b
					+s*(r+pl*(l-r))))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				s <- mixt(p)
				t <- exp(p[np])
				pl <- plaplace((y[,1]-delta/2-m)/t)
				-wt*log(s*cens*exp(-abs(y[,1]-m)/t-p[np])+
					(1-cens)*((1+s*(rc*pl-lc))*b
					+s*(r+pl*(l-r))))}
			const <- -wt*cens*log(delta/2)}},
	exponential={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- mixt(p)
				u <- exp(-(y[,1]-delta/2)/m)
				-wt*log(s*cens*(-exp(-(y[,1]+delta/2)/m)+u)
					+(1-cens)*((1+s*(rc*(1-u)-lc))*b
					+s*(r+(1-u)*(l-r))))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				s <- mixt(p)
				u <- exp(-(y[,1]-delta/2)/m)
				-wt*log(s*cens*exp(-y[,1]/m)/m
					+(1-cens)*((1+s*(rc*(1-u)-lc))*b
					+s*(r+(1-u)*(l-r))))}
			const <- -wt*cens*log(delta)}},
        gamma={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- mixt(p)
				t <- exp(p[np])
				u <- m/t
				pg <- pgamma(y[,1]-delta/2,t,u)
				-wt*log(s*cens*(pgamma(y[,1]+delta/2,t,u)-pg)
					+(1-cens)*((1+s*(rc*pg-lc))*b
					+s*(r+pg*(l-r))))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				s <- mixt(p)
				t <- exp(p[np])
				u <- m/t
				pg <- pgamma(y[,1]-delta/2,t,u)
				-wt*log(s*cens*y[,1]^(t-1)*exp(-y[,1]/u)/
					(u^t*gamma(t))
					+(1-cens)*((1+s*(rc*pg-lc))*b
					+s*(r+pg*(l-r))))}
			const <- -wt*cens*log(delta)}},
        Weibull={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- mixt(p)
				t <- exp(p[np])
				pw <- pweibull(y[,1]-delta/2,t,m)
				-wt*log(s*cens*(pweibull(y[,1]+delta/2,t,m)-pw)
					+(1-cens)*((1+s*(rc*pw-lc))*b
					+s*(r+pw*(l-r))))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				s <- mixt(p)
				t <- exp(p[np])
				v <- y[,1]/m
				u <- exp(-v^t)
				-wt*log(s*cens*t*v^(t-1)*u/m+
					(1-cens)*((1+s*(rc*(1-u)-lc))*b
					+s*(r+(1-u)*(l-r))))}
			const <- -wt*cens*log(delta)}},
        "extreme value"={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- mixt(p)
				t <- exp(p[np])
				pw <- pweibull(exp(y[,1])-delta/2,t,m)
				-wt*log(s*cens*(pweibull(exp(y[,1])+delta/2,
					t,m)-pw)+(1-cens)*((1+s*(rc*pw-lc))*b
					+s*(r+pw*(l-r))))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				s <- mixt(p)
				t <- exp(p[np])
				v <- exp(y[,1])/m
				u <- exp(-v^t)
				-wt*log(s*cens*t*v^(t-1)*u/m+
					(1-cens)*((1+s*(rc*(1-u)-lc))*b
					+s*(r+(1-u)*(l-r))))}
			const <- -wt*cens*log(delta)}},
	own={const <- 0})
fn <- function(p) sum(fcn(p))
if(fscale==1)fscale <- fn(p)
if(is.na(fscale))
	stop("Non-numerical function value: probably invalid initial values")
z0 <- nlm(fn, p=p, hessian=T, print.level=print.level, typsiz=typsiz,
	ndigit=ndigit, gradtol=gradtol, stepmax=stepmax, steptol=steptol,
	iterlim=iterlim, fscale=fscale)
z0$minimum <- z0$minimum+sum(const)
if(!is.language(lin1))cname <- paste("p",1:npl,sep="")
else {
     cname <- colnames(dm1)
     if(is.function(mu)&&length(cname)<npl)
	cname <- c(cname,paste("p",(length(cname)+1):npl,sep=""))}
if(!is.language(lin2))mname <- paste("p",1:npm,sep="")
else {
     mname <- colnames(dm2)
     if(is.function(mix))
	mname <- c(mname,paste("p",(length(mname)+1):npm,sep=""))}
fitted.values <- if(dist=="binomial"||dist=="beta binomial")
		as.vector((y[,1]+y[,2])*mu1(z0$estimate))
	else as.vector(mu1(z0$estimate))
residuals <- if(dist!="Poisson"&&dist!="negative binomial")
		y[,1]-fitted.values
	else y-fitted.values
if(npl+npm==1){
	cov <- 1/z0$hessian
	se <- sqrt(cov)}
else {
	cov <- solve(z0$hessian)
	se <- sqrt(diag(cov))}
like.comp <- as.vector(fcn(z0$estimate)+const)
if(is.function(mu))mu1 <- mu
if(is.function(mix))mixt <- mix
z1 <- list(
	call=call,
	delta=delta,
	dist=dist,
	likefn=fcn,
	mu=mu1,
	mix=mixt,
	linear=list(lin1,lin2),
	prior.weights=wt,
	censor=censor,
	maxlike=z0$minimum,
	fitted.values=fitted.values,
	residuals=residuals,
	like.comp=like.comp,
	aic=z0$minimum+np,
	df=sum(wt)-np,
	coefficients=z0$estimate,
	cname=cname,
	mname=mname,
	npl=npl,
	npm=npm,
	se=se,
	cov=cov,
	corr=cov/(se%o%se),
	gradient=z0$gradient,
	iterations=z0$iterations,
	code=z0$code)
class(z1) <- "fmr"
return(z1)}

residuals.fmr <- function(z) z$residuals
fitted.values.fmr <- function(z) z$fitted.values
coefficients.fmr <- function(z) z$coefficients
weights.fmr <- function(z) z$prior.weights
df.residual.fmr <- function(z) z$df
deviance.fmr <- function(z) 2*z$maxlike

print.fmr <- function(z) {
	sht <- z$dist!="binomial"&&z$dist!="Poisson"&&z$dist!="exponential"&&z$dist!="geometric"
	np <- z$npl+z$npm+sht
	cat("\nCall:\n",deparse(z$call),"\n\n",sep="")
	if(z$code>2)cat("Warning: no convergence - error",z$code,"\n\n")
	if(sht)cat("censored ")
	cat(z$dist,"distribution\n\n")
	t <- deparse(z$likefn)
	cat("Log likelihood function:",t[2:length(t)],"",sep="\n")
	t <- deparse(z$mu)
	cat("Location function:",t[2:length(t)],sep="\n")
	if(is.language(z$linear[[1]]))
		cat("Linear part: ",deparse(z$linear[[1]]),"\n")
	t <- deparse(z$mix)
	cat("\nMixture function:",t[2:length(t)],sep="\n")
	if(is.language(z$linear[[2]]))
		cat("Linear part: ",deparse(z$linear[[2]]),"\n")
	cat("\n-Log likelihood   ",z$maxlike,"\n")
	cat("Degrees of freedom",z$df,"\n")
	cat("AIC               ",z$aic,"\n")
	cat("Iterations        ",z$iterations,"\n\n")
	if(z$npl>0){
		cat("Location parameters:\n")
		coef.table <- cbind(z$coefficients[1:z$npl], z$se[1:z$npl])
		dimnames(coef.table) <- list(z$cname, c("estimate", "se"))
		print.default(coef.table, digits=4, print.gap=2)}
	if(z$npm>0){
		cat("\nMixture parameters:\n")
		coef.table <- cbind(z$coefficients[(z$npl+1):(np-sht)],z$se[(z$npl+1):(np-sht)])
		dimnames(coef.table) <- list(z$mname, c("estimate", "se"))
		print.default(coef.table, digits=4, print.gap=2)}
	if(sht){
		cat("\nShape parameter:\n")
		coef.table <- cbind(z$coefficients[np],z$se[np])
		dimnames(coef.table) <- list(1, c("estimate", "se"))
		print.default(coef.table, digits=4, print.gap=2)}
	if(np>1){
		cat("\nCorrelations:\n")
		dimnames(z$corr) <- list(seq(1,np),seq(1,np))
		print.default(z$corr, digits=4)}
	invisible(z)}
gnlr <- function(y, dist="normal", pmu=NULL, pshape=NULL, mu=NULL,
	shape=NULL, linear=NULL, exact=F, wt=1, delta=1, shfn=F,
	print.level=0, typsiz=abs(p), ndigit=10, gradtol=0.00001,
	stepmax=10*sqrt(p%*%p), steptol=0.00001, iterlim=100, fscale=1){
call <- sys.call()
if(!missing(dist)&&!is.function(dist)){
	dist <- match.arg(dist,c("binomial","beta binomial","double binomial",
	"mult binomial","Poisson","negative binomial","double Poisson",
	"mult Poisson","Consul","logarithmic","geometric","normal",
	"inverse Gauss","logistic","exponential","gamma","Weibull",
	"extreme value","Cauchy","Student t","Laplace"))}
if(!missing(pmu))npl <- length(pmu)
else npl <- 0
if(!missing(pshape))nps <- length(pshape)
else nps <- 0
np <- npl+nps
p <- c(pmu,pshape)
if(np<1)stop("At least one parameter must be estimated")
if(is.function(dist)){
	fcn <- dist
	dist <- "own"}
if(dist=="binomial"||dist=="double binomial"||dist=="beta binomial"||dist=="mult binomial"){
	if(length(dim(y))!=2||ncol(y)!=2)
		stop(paste("Two column matrix required for response: successes and failures"))
	if(any(y<0))stop("All response values must be positive")
	n <- nrow(y)
	nn <- y[,1]+y[,2]
	censor <- F}
else {
	if(is.list(y)&&!is.null(class(y))&&class(y)=="response"){
		if(is.null(y$censor))y <- y$y
		else y <- cbind(y$y,y$censor)}
	censor <- length(dim(y))==2&&ncol(y)==2
	if(!censor){
		if(!is.vector(y))stop("y must be a vector")
		n <- length(y)
		if(dist=="double Poisson"||dist=="mult Poisson"){
			my <- min(3*max(y),100)}}}
if((dist=="inverse Gauss"||dist=="exponential"||dist=="gamma"||
	dist=="Weibull"||dist=="extreme value")&&((censor&&any(y[,1]<=0))||
	(!censor&&any(y<=0))))stop("All response values must be > 0")
if((dist=="Poisson"||dist=="negative binomial"||
	dist=="double Poisson"||dist=="mult Poisson")&&(any(y<0)))
	stop("All response values must be >= 0")
if(dist=="logarithmic"&&any(y[wt>0]<1))
	stop("All response values must be integers > 0")
if(censor){
	n <- nrow(y)
	y[,2] <- as.integer(y[,2])
	if(any(y[,2]!=-1&y[,2]!=0&y[,2]!=1))
		stop("Censor indicator must be -1s, 0s, and 1s")
	cc <- ifelse(y[,2]==1,1,0)
	rc <- ifelse(y[,2]==0,1,ifelse(y[,2]==-1,-1,0))
	lc <- ifelse(y[,2]==-1,0,1)
	if(any(delta<=0&y[,2]==1))
		stop("All deltas for uncensored data must be positive")
	else {
		delta <- ifelse(delta<=0,0.000001,delta)
		delta <- ifelse(y[,1]-delta/2<=0,delta-0.00001,delta)}}
else {
	if(min(delta)<=0)stop("All deltas for must be positive")}
if(length(wt)==1)wt <- rep(wt,n)
else if(length(wt)!=n)stop("wt must be the same length as the other variables")
if(min(wt)<0)stop("All weights must be non-negative")
if(length(delta)==1)delta <- rep(delta,n)
else if(length(delta)!=n)stop("delta must be the same length as the other variables")
lin1 <- lin2 <- NULL
if(is.list(linear)){
	lin1 <- linear[[1]]
	lin2 <- linear[[2]]
	lsh <- is.language(lin2)}
else {
	lin1 <- linear
	lin2 <- NULL
	lsh <- FALSE}
if(is.language(mu))lin1 <- mu
if(is.language(shape)){
	lin2 <- shape
	lsh <- is.language(lin2)}
nlp <- npl
if(is.language(lin1)){
	mt <- terms(lin1)
	if(is.numeric(mt[[2]])){
		dm1 <- matrix(1)
		colnames(dm1) <- "(Intercept)"
		npt1 <- 1
		if(!is.function(mu)){
			mu1 <- function(p) p[1]*rep(1,n)
			nlp <- 1}
		else mu1 <- function(p) mu(p, p[1]*rep(1,n))}
	else {
		mf <- model.frame(mt,sys.frame(sys.parent()),na.action=na.fail)
		dm1 <- model.matrix(mt, mf)
		npt1 <- ncol(dm1)
		if(!is.function(mu)){
			mu1 <- function(p) dm1%*%p[1:npt1]
			nlp <- npt1}
		else mu1 <- function(p) mu(p, dm1%*%p[1:npt1])}
	if(npl<npt1)stop("Not enough initial estimates for mu")}
else if(!is.function(mu)){
	mu1 <- function(p) p[1]*rep(1,n)
	nlp <- 1}
else {
	mu1 <- mu
	if(length(mu1(p))==1)mu1 <- function(p) mu(p)*rep(1,n)}
if(nlp!=npl)
	stop("Number of initial estimates for mu does not correspond to model")
npl1 <- npl+1
nlp <- nps
if(lsh){
	mt <- terms(lin2)
	if(is.numeric(mt[[2]])){
		dm2 <- matrix(1)
		colnames(dm2) <- "(Intercept)"
		npt2 <- 1
		if(!is.function(shape)){
			sh1 <- function(p) p[npl1]*rep(1,n)
			nlp <- 1}
		else sh1 <- function(p) shape(p[npl1:np], p[npl1]*rep(1,n))}
	else {
		mf <- model.frame(mt,sys.frame(sys.parent()),na.action=na.fail)
		dm2 <- model.matrix(mt, mf)
		npt2 <- ncol(dm2)
		if(!is.function(shape)){
			sh1 <- function(p) dm2%*%p[npl1:(npl1+npt2-1)]
			nlp <- npt2}
		else {
			if(shfn)sh1 <- function(p) shape(p[npl1:np], dm2%*%p[npl1:(npl1+npt2-1)], mu1(p))
			else sh1 <- function(p) shape(p[npl1:np], dm2%*%p[npl1:(npl1+npt2-1)])}}
	if(nps<npt2)stop("Not enough initial estimates for shape")}
else if(!is.function(shape)&&dist!="binomial"&&dist!="Poisson"&&
	dist!="exponential"&&dist!="geometric"&&dist!="logarithmic"){
	sh1 <- function(p) p[npl1]*rep(1,n)
	nlp <- 1}
else {
	if(shfn)sh1 <- function(p) shape(p[npl1:np], mu1(p))
	else sh1 <- function(p) shape(p[npl1:np])}
if(nlp!=nps)
	stop("Number of initial estimates for shape does not correspond to model")
if(!is.numeric(mu1(pmu)))
	stop("The location function must return numerical values")
if(dist!="binomial"&&dist!="Poisson"&&dist!="exponential"&&
	dist!="geometric"&&dist!="logarithmic"&&!is.numeric(sh1(p)))
	stop("The shape function must return numerical values")
if (!censor){
	ret <- switch(dist,
	binomial={
		fcn <- function(p) {
			m <- mu1(p)
			-wt*(y[,1]*log(m)+y[,2]*log(1-m))}
		const <- -wt*lchoose(nn,y[,1])},
	"beta binomial"={
		fcn <- function(p) {
			m <- mu1(p)
			s <- exp(sh1(p))
			t <- s*m
			u <- s*(1-m)
			-wt*(lbeta(y[,1]+t,y[,2]+u)-lbeta(t,u))}
		const <- -wt*lchoose(nn,y[,1])},
	"double binomial"={
		fcn <- function(p) {
			-wt*log(.C("ddb",as.integer(y[,1]),as.integer(y[,2]),
				as.double(mu1(p)),as.double(exp(sh1(p))),
				as.integer(nrow(y)),as.integer(1),
				res=double(nrow(y)))$res)}
		const <- 0},
	"mult binomial"={
		fcn <- function(p) {
			-wt*log(.C("dmb",as.integer(y[,1]),as.integer(y[,2]),
				as.double(mu1(p)),as.double(exp(sh1(p))),
				as.integer(nrow(y)),as.integer(1),
				res=double(nrow(y)))$res)}
		const <- 0},
	Poisson={
		fcn <- function(p) {
			m <- mu1(p)
			-wt*(-m+y*log(m))}
		const <- wt*lgamma(y+1)},
	"negative binomial"={
		fcn <- function(p) {
			m <- mu1(p)
			t <- sh1(p)
			s <- exp(t)
			-wt*(lgamma(y+s)-lgamma(s)+s*t+y*log(m)
				-(y+s)*log(s+m))}
		const <- wt*lgamma(y+1)},
	"double Poisson"={
		fcn <- function(p) {
			-wt*log(.C("ddp",as.integer(y),as.integer(my),
				as.double(mu1(p)),as.double(exp(sh1(p))),
				as.integer(length(y)),as.integer(1),
				res=double(length(y)))$res)}
		const <- 0},
	"mult Poisson"={
		fcn <- function(p) {
			-wt*log(.C("dmp",as.integer(y),as.integer(my),
				as.double(mu1(p)),as.double(exp(sh1(p))),
				as.integer(length(y)),as.integer(1),
				res=double(length(y)))$res)}
		const <- 0},
	Consul={
		fcn <- function(p) {
			m <- mu1(p)
			t <- sh1(p)
			s <- exp(t)
			-wt*(log(m)-(m+y*(s-1))/s+(y-1)*log(m+y*(s-1))-y*t)}
		const <- wt*lgamma(y+1)},
	logarithmic={
		fcn <- function(p) {
			m <- exp(mu1(p))
			m <- m/(1+m)
			-wt*(y*log(m)-log(y)-log(-log(1-m)))}
		const <- 0},
	geometric={
		fcn <- function(p) {
			m <- mu1(p)
			-wt*(y*log(m)-(y+1)*log(1+m))}
		const <- 0},
	normal={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p)/2)
				-wt*log(pnorm(y+delta/2,m,s)
					-pnorm(y-delta/2,m,s))}
			const <- 0}
		else {
			fcn <- function(p) {
				t <- sh1(p)
				wt*(t+(y-mu1(p))^2/exp(t))/2}
			const <- wt*(log(2*pi)/2-log(delta))}},
        "inverse Gauss"={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p))
				-wt*log(pinvgauss(y+delta/2,m,s)-pinvgauss(y-delta/2,m,s))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				t <- sh1(p)
				wt*(t+(y-m)^2/(y*exp(t)*m^2))/2}
			const <- wt*(log(2*pi*y^3)/2-log(delta))}},
	logistic={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p))*sqrt(3)/pi
				-wt*log(plogis(y+delta/2,m,s)
					-plogis(y-delta/2,m,s))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				t <- sh1(p)
				s <- exp(t)*sqrt(3)/pi
				wt*((y-m)/s+t+2*log(1+exp(-(y-m)/s)))}
			const <- -wt*(log(pi/sqrt(3))+log(delta))}},
        Student={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p))
				-wt*log(pt(y+delta/2-m,s)-pt(y-delta/2-m,s))}
			const <- 0}
		else {
			fcn <- function(p) {
				t <- sh1(p)
				s <- exp(t)
				-wt*(lgamma((s+1)/2)-lgamma(s/2)-t/2
					-((s+1)/2)*log(1+(y-mu1(p))^2/s))}
			const <- wt*(log(pi)/2-log(delta))}},
	Cauchy={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p)/2)
				-wt*log(pcauchy(y+delta/2,m,s)
					-pcauchy(y-delta/2,m,s))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p)/2)
				wt*log(s*(1+(y-m)^2/s^2))}
			const <- -wt*log(delta/pi)}},
        Laplace={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p))
				-wt*log(plaplace((y+delta/2-m)/s)
					-plaplace((y-delta/2-m)/s))}
			const <- 0}
		else {
			fcn <- function(p) {
				t <- sh1(p)
				wt*(abs(y-mu1(p))/exp(t)+t)}
			const <- -wt*log(delta/2)}},
        exponential={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				-wt*log(-exp(-(y+delta/2)/m)
					+exp(-(y-delta/2)/m))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				wt*(log(m)+y/m)}
			const <- -wt*log(delta)}},
        gamma={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p))
				u <- m/s
				-wt*log(pgamma(y+delta/2,s,u)
					-pgamma(y-delta/2,s,u))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				t <- sh1(p)
				s <- exp(t)
				-wt*(s*(t-log(m)-y/m)+(s-1)*log(y)-lgamma(s))}
			const <- -wt*log(delta)}},
        Weibull={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p))
				-wt*log(pweibull(y+delta/2,s,m)
					-pweibull(y-delta/2,s,m))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				t <- sh1(p)
				s <- exp(t)
				-wt*(t+(s-1)*log(y)-s*log(m)-(y/m)^s)}
			const <- -wt*log(delta)}},
        "extreme value"={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p))
				ey <- exp(y[,1])
				-wt*log(pweibull(ey+ey*delta/2,s,m)
					-pweibull(ey-ey*delta/2,s,m))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				t <- sh1(p)
				s <- exp(t)
				-wt*(t+s*y-s*log(m)-(exp(y)/m)^s)}
			const <- -wt*log(delta)}},
	own={ const <- 0})}
else {
	ret <- switch(dist,
	Poisson={
		fcn <- function(p) {
			m <- mu1(p)
			-wt*(cc*(-m+y[,1]*log(m))+
				log(lc-rc*ppois(y[,1],m)))}
		const <- wt*cc*lgamma(y[,1]+1)},
	"negative binomial"={
		fcn <- function(p) {
			m <- mu1(p)
			t <- sh1(p)
			s <- exp(t)
			-wt*(cc*(lgamma(y[,1]+s)-lgamma(s)
				+s*t+y[,1]*log(m)-(y[,1]+s)*log(s+m))+
				log(lc-rc*pnbinom(y[,1],s,1/(1+m/s))))}
		const <- wt*cc*lgamma(y[,1]+1)},
	geometric={
		fcn <- function(p) {
			m <- mu1(p)
			-wt*(cc*(y[,1]*log(m)-(y[,1]+1)*log(1+m))+
				log(lc-rc*pgeom(y[,1],1/(1+m))))}
		const <- 0},
	normal={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p)/2)
				pn <- pnorm(y[,1]-delta/2,m,s)
				-wt*(cc*log(pnorm(y[,1]+delta/2,m,s)-pn)
					+log(lc-rc*pn))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				t <- sh1(p)
				s <- exp(t)
				-wt*(cc*(-(t+(y[,1]-m)^2/s)/2)+log(lc-rc
					*pnorm(y[,1]-delta/2,m,sqrt(s))))}
			const <- wt*cc*(log(2*pi)/2-log(delta))}},
        "inverse Gauss"={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p))
				pit <- pinvgauss(y[,1]-delta/2,m,s)
				-wt*(cc*log(pinvgauss(y[,1]+delta/2,m,s)-pit)
					+log(lc-rc*pit))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				t <- sh1(p)
				s <- exp(t)
				-wt*(cc*(-(t+(y[,1]-m)^2/(y[,1]*s*m^2))/2)+
					log(lc-rc*pinvgauss(y[,1]-delta/2,m,s)))}
			const <- wt*cc*(log(2*pi*y[,1]^3)/2-log(delta))}},
	logistic={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p))*sqrt(3)/pi
				pl <- plogis(y[,1]-delta/2,m,s)
				-wt*(cc*log(plogis(y[,1]+delta/2,m,s)-pl)
					+log(lc-rc*pl))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p))*sqrt(3)/pi
				y1 <- (y[,1]-m)/s
				-wt*(cc*(-y1-log(s)-2*log(1+exp(-y1)))
					+log(lc-rc*plogis(y[,1]-delta/2,m,s)))}
			const <- -wt*cc*log(delta)}},
        Student={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p))
				ps <- pt(y[,1]-delta/2-m,s)
				-wt*(cc*log(pt(y[,1]+delta/2-m,s)-ps)
					+log(lc-rc*ps))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				t <- sh1(p)
				s <- exp(t)
				-wt*(cc*(lgamma((s+1)/2)-lgamma(s/2)-t/2
					-((s+1)/2)*log(1+(y[,1]-m)^2/s))
					+log(lc-rc*pt(y[,1]-delta/2-m,s)))}
			const <- wt*cc*(log(pi)/2-log(delta))}},
	Cauchy={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p)/2)
				pc <- pcauchy(y[,1]-delta/2,m,s)
				-wt*(cc*log(pcauchy(y[,1]+delta/2,m,s)-pc)
					+log(lc-rc*pc))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p)/2)
				-wt*(-cc*log(s*(1+(y[,1]-m)^2/s^2))
					+log(lc-rc*pcauchy(y[,1]-delta/2
					,m,s)))}
			const <- -wt*cc*log(delta/pi)}},
        Laplace={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p))
				pl <- plaplace((y[,1]-delta/2-m)/s)
				-wt*(cc*log(plaplace((y[,1]+delta/2-m)/s)-pl)
					+log(lc-rc*pl))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				t <- sh1(p)
				s <- exp(t)
				-wt*(cc*(-abs(y[,1]-m)/s-t)+log(lc-rc
					*plaplace((y[,1]-delta/2-m)/s)))}
			const <- -wt*cc*log(delta/2)}},
	exponential={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				pe <- exp(-(y[,1]-delta/2)/m)
				-wt*(cc*log(-exp(-(y[,1]+delta/2)/m)
					+pe)+log(lc-rc*(1-pe)))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				-wt*(cc*(-log(m)-y[,1]/m)+log(lc-rc*
					(1-exp(-(y[,1]-delta/2)/m))))}
			const <- -wt*cc*log(delta)}},
        gamma={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p))
				u <- m/s
				pg <- pgamma(y[,1]-delta/2,s,u)
				-wt*(cc*log(pgamma(y[,1]+delta/2,s,u)-pg)
					+log(lc-rc*pg))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				t <- sh1(p)
				s <- exp(t)
				-wt*(cc*(s*(t-log(m)-y[,1]/m)+(s-1)*log(y[,1])
					-lgamma(s))+log(lc-rc
					*pgamma(y[,1]-delta/2,s,m/s)))}
			const <- -wt*cc*log(delta)}},
        Weibull={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p))
				pw <- pweibull(y[,1]-delta/2,s,m)
				-wt*(cc*log(pweibull(y[,1]+delta/2,s,m)-pw)
					+log(lc-rc*pw))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				t <- sh1(p)
				s <- exp(t)
				-wt*(cc*(t+(s-1)*log(y[,1])-s*log(m)
					-(y[,1]/m)^s)+log(lc-rc*
					pweibull(y[,1]-delta/2,s,m)))}
			const <- -wt*cc*log(delta)}},
        "extreme value"={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p))
				ey <- exp(y[,1])
				pw <- pweibull(ey-ey*delta/2,s,m)
				-wt*(cc*log(pweibull(ey+ey*delta/2,s,m)-pw)
					+log(lc-rc*pw))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				t <- sh1(p)
				s <- exp(t)
				ey <- exp(y[,1])
				-wt*(cc*(t+s*y[,1]-s*log(m)-(ey/m)^s)+log(lc-
					rc*pweibull(ey-ey*delta/2,s,m)))}
			const <- -wt*cc*log(delta)}},
	own={const <- 0})}
fn <- function(p) sum(fcn(p))
if(fscale==1)fscale <- fn(p)
if(is.na(fscale))
	stop("Non-numerical function value: probably invalid initial values")
z0 <- nlm(fn, p=p, hessian=T, print.level=print.level, typsiz=typsiz,
	ndigit=ndigit, gradtol=gradtol, stepmax=stepmax, steptol=steptol,
	iterlim=iterlim, fscale=fscale)
z0$minimum <- z0$minimum+sum(const)
if(!is.language(lin1))cname <- paste("p",1:npl,sep="")
else {
     cname <- colnames(dm1)
     if(is.function(mu)&&length(cname)<npl)
	cname <- c(cname,paste("p",(length(cname)+1):npl,sep=""))}
if(!is.language(lin2))sname <- paste("p",1:nps,sep="")
else {
     sname <- colnames(dm2)
     if(is.function(shape))
	sname <- c(sname,paste("p",(length(sname)+1):nps,sep=""))}
fitted.values <- if(dist=="binomial"||dist=="beta binomial"||dist=="double binomial"||dist=="mult binomial")
		as.vector((y[,1]+y[,2])*mu1(z0$estimate))
	else as.vector(mu1(z0$estimate))
residuals <- if(dist=="binomial"||dist=="beta binomial"||dist=="double binomial"||dist=="mult binomial"||censor)
		y[,1]-fitted.values
	else y-fitted.values
if(npl+nps==1){
	cov <- 1/z0$hessian
	se <- sqrt(cov)}
else {
	a <- qr(z0$hessian)
	if(a$rank==npl+nps)cov <- solve(z0$hessian)
	else cov <- matrix(NA,ncol=npl+nps,nrow=npl+nps)
	se <- sqrt(diag(cov))}
like.comp <- as.vector(fcn(z0$estimate)+const)
if(is.function(mu))mu1 <- mu
if(is.function(shape))sh1 <- shape
z1 <- list(
	call=call,
	delta=delta,
	dist=dist,
	likefn=fcn,
	mu=mu1,
	shape=sh1,
	linear=list(lin1,lin2),
	prior.weights=wt,
	censor=censor,
	maxlike=z0$minimum,
	fitted.values=fitted.values,
	residuals=residuals,
	like.comp=like.comp,
	aic=z0$minimum+np,
	df=sum(wt)-np,
	coefficients=z0$estimate,
	cname=cname,
	sname=sname,
	npl=npl,
	nps=nps,
	se=se,
	cov=cov,
	corr=cov/(se%o%se),
	gradient=z0$gradient,
	iterations=z0$iterations,
	code=z0$code)
class(z1) <- "gnlr"
return(z1)}

residuals.gnlr <- function(z) z$residuals
fitted.values.gnlr <- function(z) z$fitted.values
coefficients.gnlr <- function(z) z$coefficients
weights.gnlr <- function(z) z$prior.weights
df.residual.gnlr <- function(z) z$df
deviance.gnlr <- function(z) 2*z$maxlike

print.gnlr <- function(z) {
	np <- z$npl+z$nps
	cat("\nCall:\n",deparse(z$call),"\n\n",sep="")
	if(z$code>2)cat("Warning: no convergence - error",z$code,"\n\n")
	if(z$censor)cat("censored ")
	cat(z$dist,"distribution\n\n")
	t <- deparse(z$likefn)
	cat("Log likelihood function:",t[2:length(t)],"",sep="\n")
	t <- deparse(z$mu)
	cat("Location function:",t[2:length(t)],sep="\n")
	if(is.language(z$linear[[1]]))
		cat("Linear part: ",deparse(z$linear[[1]]),"\n")
	if(z$dist!="binomial"&&z$dist!="Poisson"&&z$dist!="geometric"&&
		z$dist!="exponential"){
		t <- deparse(z$shape)
		cat("\nLog shape function:",t[2:length(t)],sep="\n")
		if(!is.null(z$linear[[2]])&&is.language(z$linear[[2]]))
			cat("Linear part: ",deparse(z$linear[[2]]),"\n")}
	cat("\n-Log likelihood   ",z$maxlike,"\n")
	cat("Degrees of freedom",z$df,"\n")
	cat("AIC               ",z$aic,"\n")
	cat("Iterations        ",z$iterations,"\n\n")
	if(z$npl>0){
		cat("Location parameters:\n")
		coef.table <- cbind(z$coefficients[1:z$npl],z$se[1:z$npl])
		dimnames(coef.table) <- list(z$cname, c("estimate", "se"))
		print.default(coef.table, digits=4, print.gap=2)}
	if(z$nps>0){
		cat("\nShape parameters:\n")
		coef.table <- cbind(z$coefficients[(z$npl+1):np],z$se[(z$npl+1):np])
		dimnames(coef.table) <- list(z$sname, c("estimate", "se"))
		print.default(coef.table, digits=4, print.gap=2)}
	if(np>1){
		cat("\nCorrelations:\n")
		dimnames(z$corr) <- list(seq(1,np),seq(1,np))
		print.default(z$corr, digits=4)}
	invisible(z)}

wr <- function(formula)
{
	mt <- terms(formula)
	data <- sys.frame(sys.parent())
	mf <- model.frame(mt, data, na.action=na.fail)
	x <- model.matrix(mt, mf)
	y <- model.response(mf, "numeric")
	z <- list(response=y, design=x)
	z
}
gnlr3 <- function(y, dist="normal", pmu=NULL, pshape=NULL,
	pfamily=NULL, mu=NULL, shape=NULL, family=NULL,
	linear=NULL, exact=F, wt=1, delta=1, print.level=0,
	typsiz=abs(p), ndigit=10, gradtol=0.00001, stepmax=10*sqrt(p%*%p),
	steptol=0.00001, iterlim=100, fscale=1){
pweib <- function(y,s,m,f) 1-(1-exp(-(y/m)^s))^f
plog <- function(y,m,s,f) (1+exp(-(y-m)/s))^-f
phjorth <- function(y,m,s,f) 1-(1+s*y)^(-f/s)*exp(-(y/m)^2/2)
pburr <- function(y,m,s,f) 1-(1+(y/m)^s)^-f
#pinvg <- function(y,m,s,f){
#	ll <- function(y) m^-f*y^(f-1)*exp(-s*(m/y+y/m)/2)/(2*bessel(s,f)$rk)
#	z <- rep(0,length(y))
#	for(i in 1:length(y))z[i] <- int(ll,"i",y[i])
#	z}
call <- sys.call()
if(!missing(dist)&&!is.function(dist)){
	dist <- match.arg(dist,c("normal","inverse Gauss","logistic",
	"Hjorth","gamma","Burr","Weibull","extreme value","elliptic"))}
if(!missing(pmu))npl <- length(pmu)
else npl <- 0
if(!missing(pshape))nps <- length(pshape)
else nps <- 0
if(!missing(pfamily))npf <- length(pfamily)
else npf <- 0
np <- npl+nps+npf
if(np<1)stop("At least one parameter must be estimated")
if(is.function(dist)){
	fcn <- dist
	dist <- "own"}
if(is.list(y)&&!is.null(class(y))&&class(y)=="response"){
	if(is.null(y$censor))y <- y$y
	else y <- cbind(y$y,y$censor)}
censor <- length(dim(y))==2&&ncol(y)==2
if(censor){
	n <- nrow(y)
	y[,2] <- as.integer(y[,2])
	if(any(y[,2]!=-1&y[,2]!=0&y[,2]!=1))
		stop("Censor indicator must be -1s, 0s, and 1s")
	cc <- ifelse(y[,2]==1,1,0)
	rc <- ifelse(y[,2]==0,1,ifelse(y[,2]==-1,-1,0))
	lc <- ifelse(y[,2]==-1,0,1)
	if(delta<=0&y[,2]==1)
		stop("All deltas for uncensored data must be positive")
	else {
		delta <- ifelse(delta<=0,0.000001,delta)
		delta <- ifelse(y[,1]-delta/2<=0,delta-0.00001,delta)}}
else {
	if(!is.vector(y))stop("y must be a vector")
	n <- length(y)
	if(min(delta)<=0)stop("All deltas for must be positive")}
if((dist!="logistic"&&dist!="elliptic")&&((censor&&any(y[,1]<=0))||
	(!censor&&any(y<=0))))stop("All response values must be > 0")
if(min(wt)<0)stop("All weights must be non-negative")
if(length(wt)==1)wt <- rep(wt,n)
if(length(delta)==1)delta <- rep(delta,n)
lin <- list(NULL,NULL,NULL)
if(!is.null(linear)){
	if(is.list(linear))lin <- linear
	else lin[[1]] <- linear}
if(is.language(mu))lin[[1]] <- mu
if(is.language(shape))lin[[2]] <- shape
if(is.language(family))lin[[3]] <- family
nlp <- npl
if(is.language(lin[[1]])){
	mt <- terms(lin[[1]])
	if(is.numeric(mt[[2]])){
		dm1 <- matrix(1)
		colnames(dm1) <- "(Intercept)"
		npt1 <- 1
		if(!is.function(mu)){
			mu1 <- function(p) p[1]*rep(1,n)
			nlp <- 1}
		else mu1 <- function (p) mu(p, p[1]*rep(1,n))}
	else {
		mf <- model.frame(mt,sys.frame(sys.parent()),na.action=na.fail)
		dm1 <- model.matrix(mt, mf)
		npt1 <- ncol(dm1)
		if(!is.function(mu)){
			mu1 <- function(p) dm1 %*% p[1:npt1]
			nlp <- npt1}
		else mu1 <- function (p) mu(p, dm1 %*% p[1:npt1])}
	if(npl<npt1)stop("Not enough initial estimates for mu")}
else if(!is.function(mu)){
	mu1 <- function(p) p[1]*rep(1,n)
	nlp <- 1}
else mu1 <- mu
if(nlp!=npl)
	stop("Number of initial estimates for mu does not correspond to model")
nlp <- nps
if(is.language(lin[[2]])){
	mt <- terms(lin[[2]])
	if(is.numeric(mt[[2]])){
		dm2 <- matrix(1)
		colnames(dm2) <- "(Intercept)"
		npt2 <- 1
		if(!is.function(shape)){
			sh1 <- function(p) p[1]*rep(1,n)
			nlp <- 1}
		else sh1 <- function(p) shape(p, p[1]*rep(1,n))}
	else {
		mf <- model.frame(mt,sys.frame(sys.parent()),na.action=na.fail)
		dm2 <- model.matrix(mt, mf)
		npt2 <- ncol(dm2)
		if(!is.function(shape)){
			sh1 <- function(p)dm2 %*% p[1:npt2]
			nlp <- npt2}
		else sh1 <- function(p)shape(p, dm2 %*% p[1:npt2])}
	if(nps<npt2)stop("Not enough initial estimates for shape")}
else if(!is.function(shape)){
	sh1 <- function(p) p[1]*rep(1,n)
	nlp <- 1}
else sh1 <- shape
if(nlp!=nps)
	stop("Number of initial estimates for shape does not correspond to model")
nlp <- npf
if(is.language(lin[[3]])){
	mt <- terms(lin[[3]])
	if(is.numeric(mt[[2]])){
		dm3 <- matrix(1)
		colnames(dm3) <- "(Intercept)"
		npt3 <- 1
		if(!is.function(family)){
			fa1 <- function(p) p[1]*rep(1,n)
			nlp <- 1}
		else fa1 <- function(p) family(p,p[1]*rep(1,n))}
	else {
		mf <- model.frame(mt,sys.frame(sys.parent()),na.action=na.fail)
		dm3 <- model.matrix(mt, mf)
		npt3 <- ncol(dm3)
		if(!is.function(family)){
			fa1 <- function(p) dm3 %*% p[1:npt3]
			nlp <- npt3}
		else fa1 <- function(p)family(p, dm3 %*% p[1:npt3])}
	if(npf<npt3)stop("Not enough initial estimates for family")}
else if(!is.function(family)){
	fa1 <- function(p) p[1]*rep(1,n)
	nlp <- 1}
else fa1 <- family
if(nlp!=npf)
	stop("Number of initial estimates for family does not correspond to model")
if(!is.numeric(mu1(pmu)))stop("The mean function must return numerical values")
if(!is.numeric(sh1(pshape)))
	stop("The shape function must return numerical values")
if(!is.numeric(fa1(pfamily)))
	stop("The family function must return numerical values")
p <- c(pmu,pshape,pfamily)
npl1 <- npl+1
np1 <- npl+nps
nps1 <- np1+1
if (!censor){
	ret <- switch(dist,
	normal={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p[npl1:np1])/2)
				f <- fa1(p[nps1:np])
				-wt*log(pnorm((y+delta/2)^f,m,s)
					-pnorm((y-delta/2)^f,m,s))}
			const <- 0}
		else {
			fcn <- function(p) {
				t <- sh1(p[npl1:np1])
				f <- fa1(p[nps1:np])
				wt*((t+(y^f/f-mu1(p))^2/exp(t))/2
					-(f-1)*log(y))}
			const <- wt*(log(2*pi)/2-log(delta))}},
	elliptic={
		if(exact){stop("Not yet available")}
		else {
			fcn <- function(p) {
				t <- 0.5*sh1(p[npl1:np1])
				f <- exp(fa1(p[nps1:np]))
				b <- 1+1/(2*f)
				wt*(t+(abs(y-mu1(p))/exp(t))^(2*f)/2+
					lgamma(b)+b*log(2))}
			const <- -wt*log(delta)}},
	"inverse Gauss"={
		if(exact)stop("Not yet available")
		else {
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p[npl1:np1]))
				f <- fa1(p[nps1:np])
				b <- log(2*bessel(1/(s*m),f)$rk)
				-wt*(-f*log(m)+(f-1)*log(y)-b-(1/y+y/m^2)/(2*s))}
			const <- -wt*log(delta)}},
	logistic={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p[npl1:np1]))*sqrt(3)/pi
				f <- exp(fa1(p[nps1:np]))
				-wt*log(plog(y+delta/2,m,s,f)
					-plog(y-delta/2,m,s,f))}
			const <- 0}
		else {
			fcn <- function(p) {
				t <- sh1(p[npl1:np1])
				m <- (y-mu1(p))/exp(t)*sqrt(3)/pi
				wt*(-fa1(p[nps1:np])+m+t+(exp(fa1(p[nps1:np]))
					+1)*log(1+exp(-m)))}
			const <- -wt*(log(pi/sqrt(3))+log(delta))}},
	Hjorth={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p[npl1:np1]))
				f <- fa1(p[nps1:np])
				-wt*log(phjorth(y+delta/2,m,s,f)-
					phjorth(y-delta/2,m,s,f))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p[npl1:np1]))
				f <- fa1(p[nps1:np])
				-wt*(-f*log(1+s*y)/s-(y/m)^2/2+
					log(y/m^2+f/(1+s*y)))}
			const <- -wt*log(delta)}},
        gamma={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p[npl1:np1]))
				f <- exp(fa1(p[nps1:np]))
				u <- (m/s)^f
				-wt*log(pgamma((y+delta/2)^f,s,u)
					-pgamma((y-delta/2)^f,s,u))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				t <- sh1(p[npl1:np1])
				s <- exp(t)
				u <- fa1(p[nps1:np])
				f <- exp(u)
				v <- s*f
				-wt*(v*(t-log(m))-(s*y/m)^f+u+(v-1)*log(y)
					-lgamma(s))}
			const <- -wt*log(delta)}},
	Burr={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p[npl1:np1]))
				f <- exp(fa1(p[nps1:np]))
				-wt*log(pburr(y+delta/2,m,s,f)-
					pburr(y-delta/2,m,s,f))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p[npl1:np1]))
				f <- exp(fa1(p[nps1:np]))
				y1 <- y/m
				-wt*(log(f*s/m)+(s-1)*log(y1)
					-(f+1)*log(1+y1^s))}
			const <- -wt*log(delta)}},
        Weibull={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p[npl1:np1]))
				f <- exp(fa1(p[nps1:np]))
				-wt*log(pweib(y+delta/2,s,m,f)
					-pweib(y-delta/2,s,m,f))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				t <- sh1(p[npl1:np1])
				s <- exp(t)
				u <- fa1(p[nps1:np])
				f <- exp(u)
				y1 <- (y/m)^s
				-wt*(t+u+(s-1)*log(y)-s*log(m)+
					(f-1)*log(1-exp(-y1))-y1)}
			const <- -wt*log(delta)}},
        "extreme value"={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p[npl1:np1]))
				f <- fa1(p[nps1:np])
				y1 <- y^f/f
				ey <- exp(y1)
				jey <- y^(f-1)*ey
				-wt*log(pweibull(ey+jey*delta/2,s,m)
					-pweibull(ey-jey*delta/2,s,m))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				t <- sh1(p[npl1:np1])
				s <- exp(t)
				f <- fa1(p[nps1:np])
				y1 <- y^f/f
				-wt*(t+s*(y1-log(m))-(exp(y1)/m)^s
					+(f-1)*log(y))}
			const <- -wt*log(delta)}},
	own={ const <- 0})}
else {
	ret <- switch(dist,
	normal={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p[npl1:np1])/2)
				f <- fa1(p[nps1:np])
				pn <- pnorm((y[,1]-delta/2)^f/f,m,s)
				-wt*(cc*log(pnorm((y[,1]+delta/2)^f/f,m,s)-pn)
					+log(lc-rc*pn))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				t <- sh1(p[npl1:np1])
				s <- exp(t)
				f <- fa1(p[nps1:np])
				-wt*(cc*(-(t+(y[,1]^f/f-m)^2/s)/2+(f-1)*
					log(y[,1]))+log(lc-rc
					*pnorm((y[,1]-delta/2)^f/f,m,sqrt(s))))}
			const <- wt*cc*(log(2*pi)/2-log(delta))}},
	elliptic={
		stop("Not yet available")},
	"inverse Gauss"={
		stop("Not yet available")},
	logistic={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p[npl1:np1]))*sqrt(3)/pi
				f <- exp(fa1(p[nps1:np]))
				pl <- plog(y[,1]-delta/2,m,s,f)
				-wt*(cc*log(plog(y[,1]+delta/2,m,s,f)-pl)
					+log(lc-rc*pl))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p[npl1:np1]))*sqrt(3)/pi
				y1 <- (y[,1]-m)/s
				u <- fa1(p[nps1:np])
				f <- exp(u)
				-wt*(cc*(u-y1-log(s)-(f+1)*log(1+exp(-y1)))
					+log(lc-rc*plog(y[,1]-delta/2,m,s,f)))}
			const <- -wt*cc*log(delta)}},
	Hjorth={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p[npl1:np1]))
				f <- fa1(p[nps1:np])
				ph <- phjorth(y[,1]-delta/2,m,s,f)
				-wt*(cc*log(phjorth(y[,1]+delta/2,m,s,f)-ph)
					+log(lc-rc*ph))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p[npl1:np1]))
				f <- fa1(p[nps1:np])
				-wt*(cc*(-f*log(1+s*y[,1])/s-(y[,1]/m)^2/2+
					log(y[,1]/m^2+f/(1+s*y[,1])))+
					log(lc-rc*phjorth(y[,1]
					-delta/2,m,s,f)))}
			const <- -wt*cc*log(delta)}},
        gamma={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p[npl1:np1]))
				f <- exp(fa1(p[nps1:np]))
				u <- (m/s)^f
				pg <- pgamma((y[,1]-delta/2)^f,s,u)
				-wt*(cc*log(pgamma((y[,1]+delta/2)^f,s,u)-pg)
					+log(lc-rc*pg))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				t <- sh1(p[npl1:np1])
				s <- exp(t)
				u <- fa1(p[nps1:np])
				f <- exp(u)
				v <- s*f
				-wt*(cc*(v*(t-log(m))-(s*y[,1]/m)^f+u+(v-1)*
					log(y[,1])-lgamma(s))+log(lc-rc
					*pgamma((y[,1]-delta/2)^f,s,(m/s)^f)))}
			const <- -wt*cc*log(delta)}},
	Burr={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p[npl1:np1]))
				f <- exp(fa1(p[nps1:np]))
				pb <- pburr(y[,1]-delta/2,m,s,f)
				-wt*(cc*log(pburr(y[,1]+delta/2,m,s,f)-pb)
					+log(lc-rc*pb))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p[npl1:np1]))
				f <- exp(fa1(p[nps1:np]))
				y1 <- y[,1]/m
				-wt*(cc*(log(f*s/m)+(s-1)*log(y1)
					-(f+1)*log(1+(y1)^s))+
					log(lc-rc*pburr(y[,1]
					-delta/2,m,s,f)))}
			const <- -wt*cc*log(delta)}},
        Weibull={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p[npl1:np1]))
				f <- exp(fa1(p[nps1:np]))
				pw <- pweib(y[,1]-delta/2,s,m,f)
				-wt*(cc*log(pweib(y[,1]+delta/2,s,m,f)-pw)
					+log(lc-rc*pw))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				t <- sh1(p[npl1:np1])
				s <- exp(t)
				u <- fa1(p[nps1:np])
				f <- exp(u)
				y1 <- (y[,1]/m)^s
				-wt*(cc*(t+u+(s-1)*log(y[,1])-s*log(m)+(f-1)
					*log(1-exp(-y1))-y1)+log(lc-rc*
					pweib(y[,1]-delta/2,s,m,f)))}
			const <- -wt*cc*log(delta)}},
        "extreme value"={
		if(exact){
			fcn <- function(p) {
				m <- mu1(p)
				s <- exp(sh1(p[npl1:np1]))
				f <- fa1(p[nps1:np])
				y1 <- y[,1]^f/f
				ey <- exp(y1)
				jey <- y[,1]^(f-1)*ey
				pw <- pweibull(ey-jey*delta/2,s,m)
				-wt*(cc*log(pweibull(ey+jey*delta/2,s,m)-pw)
					+log(lc-rc*pw))}
			const <- 0}
		else {
			fcn <- function(p) {
				m <- mu1(p)
				t <- sh1(p[npl1:np1])
				s <- exp(t)
				f <- fa1(p[nps1:np])
				y1 <- y[,1]^f/f
				ey <- exp(y1)
				-wt*(cc*(t+s*(y1-log(m))-(ey/m)^s
					+(f-1)*log(y[,1]))+log(lc-rc*
					pweibull(ey-y[,1]^(f-1)*ey
					*delta/2,s,m)))}
			const <- -wt*cc*log(delta)}},
	own={const <- 0})}
fn <- function(p) sum(fcn(p))
if(fscale==1)fscale <- fn(p)
if(is.na(fscale))
	stop("Non-numerical function value: probably invalid initial values")
z0 <- nlm(fn, p=p, hessian=T, print.level=print.level, typsiz=typsiz,
	ndigit=ndigit, gradtol=gradtol, stepmax=stepmax, steptol=steptol,
	iterlim=iterlim, fscale=fscale)
z0$minimum <- z0$minimum+sum(const)
if(!is.language(lin[[1]]))cname <- paste("p",1:npl,sep="")
else {
     cname <- colnames(dm1)
     if(is.function(mu)&&length(cname)<npl)
	cname <- c(cname,paste("p",(length(cname)+1):npl,sep=""))}
if(!is.language(lin[[2]]))sname <- paste("p",1:nps,sep="")
else {
     sname <- colnames(dm2)
     if(is.function(shape))
	sname <- c(sname,paste("p",(length(sname)+1):nps,sep=""))}
if(!is.language(lin[[3]]))fname <- paste("p",1:npf,sep="")
else {
     fname <- colnames(dm3)
     if(is.function(family))
	fname <- c(fname,paste("p",(length(fname)+1):npf,sep=""))}
fitted.values <- as.vector(mu1(z0$estimate))
residuals <- y-fitted.values
if(np==1){
	cov <- 1/z0$hessian
	se <- sqrt(cov)}
else {
	cov <- solve(z0$hessian)
	se <- sqrt(diag(cov))}
like.comp <- as.vector(fcn(z0$estimate)+const)
if(is.function(mu))mu1 <- mu
if(is.function(shape))sh1 <- shape
if(is.function(family))fa1 <- family
z1 <- list(
	call=call,
	delta=delta,
	dist=dist,
	likefn=fcn,
	mu=mu1,
	shape=sh1,
	family=fa1,
	linear=lin,
	prior.weights=wt,
	censor=censor,
	maxlike=z0$minimum,
	fitted.values=fitted.values,
	residuals=residuals,
	like.comp=like.comp,
	aic=z0$minimum+np,
	df=sum(wt)-np,
	coefficients=z0$estimate,
	cname=cname,
	sname=sname,
	fname=fname,
	npl=npl,
	nps=nps,
	npf=npf,
	se=se,
	cov=cov,
	corr=cov/(se%o%se),
	gradient=z0$gradient,
	iterations=z0$iterations,
	code=z0$code)
class(z1) <- "gnlr3"
return(z1)}

residuals.gnlr3 <- function(z) z$residuals
fitted.values.gnlr3 <- function(z) z$fitted.values
coefficients.gnlr3 <- function(z) z$coefficients
weights.gnlr3 <- function(z) z$prior.weights
df.residual.gnlr3 <- function(z) z$df
deviance.gnlr3 <- function(z) 2*z$maxlike

print.gnlr3 <- function(z) {
	np1 <- z$npl+1
	np2 <- z$npl+z$nps
	np3 <- np2+1
	np <- z$npl+z$nps+z$npf
	cat("\nCall:\n",deparse(z$call),"\n\n",sep="")
	if(z$code>2)cat("Warning: no convergence - error",z$code,"\n\n")
	if(z$censor)cat("censored ")
	cat(z$dist,"distribution\n\n")
	t <- deparse(z$likefn)
	cat("Log likelihood function:",t[2:length(t)],"",sep="\n")
	t <- deparse(z$mu)
	cat("Mean function:",t[2:length(t)],sep="\n")
	if(is.language(z$linear[[1]]))
		cat("Linear part: ",deparse(z$linear[[1]]),"\n")
	t <- deparse(z$shape)
	cat("\nLog shape function:",t[2:length(t)],sep="\n")
	if(is.language(z$linear[[2]]))
		cat("Linear part: ",deparse(z$linear[[2]]),"\n")
	t <- deparse(z$family)
	cat("\n(Log) family function:",t[2:length(t)],sep="\n")
	if(is.language(z$linear[[3]]))
		cat("Linear part: ",deparse(z$linear[[3]]),"\n")
	cat("\n-Log likelihood   ",z$maxlike,"\n")
	cat("Degrees of freedom",z$df,"\n")
	cat("AIC               ",z$aic,"\n")
	cat("Iterations        ",z$iterations,"\n\n")
	if(z$npl>0){
		cat("Location parameters:\n")
		coef.table <- cbind(z$coefficients[1:z$npl], z$se[1:z$npl])
		dimnames(coef.table) <- list(z$cname, c("estimate", "se"))
		print.default(coef.table, digits=4, print.gap=2)}
	if(z$nps>0){
		cat("\nShape parameters:\n")
		coef.table <- cbind(z$coefficients[np1:np2],z$se[np1:np2])
		dimnames(coef.table) <- list(z$sname, c("estimate", "se"))
		print.default(coef.table, digits=4, print.gap=2)}
	if(z$npf>0){
		cat("\nFamily parameters:\n")
		coef.table <- cbind(z$coefficients[np3:np],z$se[np3:np])
		dimnames(coef.table) <- list(z$fname, c("estimate", "se"))
		print.default(coef.table, digits=4, print.gap=2)}
	if(np>1){
		cat("\nCorrelations:\n")
		dimnames(z$corr) <- list(seq(1,np),seq(1,np))
		print.default(z$corr, digits=4)}
	invisible(z)}
library.dynam("gnlm.so")

int <- function(f, a, b, eps=1.0e-6, k=5)
{
	ff <- function(x) f(1/x)/(x*x)
	if(!is.numeric(b)){
		if(!is.numeric(a)) z <- int1(ff, -1, 0, eps, k) +
			int1(f, -1, 1, eps, k) + int1(ff, 0, 1, eps, k)
		else {
			if(a>0)	z <- int1(ff, 0, 1/a, eps, k)
			else z <- int1(f, a, 1, eps, k) +
				int1(ff, 0, 1, eps, k)}}
	else if(!is.numeric(a)){
		if(b<0)	z <- int1(ff, 1/b, 0, eps, k)
		else z <- int1(f, -1, b, eps, k) + int1(ff, -1, 0, eps, k)}
	else z <- int1(f, a, b, eps, k)
	z
}
int1 <- function(ff, aa, bb, eps, k)
{
	z <- .C("qromo",
		ff,
		aa,
		bb,
		eps=eps,
		k=as.integer(k),
		err=integer(1),
		res=double(1))
	if(z$err>0)print("Int: non-zero error code")
	z$res
}
nlr <- function(y, mu=NULL, p=NULL, dist="normal", wt=1, delta=1,
	print.level=0, typsiz=abs(p), ndigit=10, gradtol=0.00001,
	stepmax=10*sqrt(p%*%p), steptol=0.00001, iterlim=100, fscale=1){
	call <- sys.call()
	if(!missing(dist))dist <- match.arg(dist,c("normal","inverse Gauss","gamma"))
	if(missing(p))stop("Initial parameter estimates must be supplied")
	else np <- length(p)
	if(missing(mu)||!is.function(mu))stop("A mean function must be supplied")
	switch(dist,
	normal=fn <- function(p) sum(wt*(y-mu(p))^2),
	gamma=fn <- function(p) -sum(wt*(log(y/mu(p))-(y-mu(p))/mu(p))),
	"inverse Gauss"=fn <- function(p) sum(wt*((y-mu(p))^2)/(y*mu(p)^2)))
	if(fscale==1)fscale <- fn(p)
	if(is.na(fscale))
		stop("Non-numerical function value: probably invalid initial values")
	z0 <- nlm(fn, p=p, hessian=T, print.level=print.level, typsiz=typsiz,
		ndigit=ndigit, gradtol=gradtol, stepmax=stepmax,
		steptol=steptol, iterlim=iterlim, fscale=fscale)
	if(length(delta)==1)delta <- rep(delta,length(y))
	if(length(wt)==1)wt <- rep(wt,length(y))
	n <- sum(wt)
	disp <- z0$minimum/n
	p <- z0$estimate
	switch(dist,
	normal=maxlike <- length(y)*(log(2*pi*disp)+1)/2,
	gamma=maxlike <- (sum(wt*(y/mu(p)+log(mu(p))-log(y)))+n*log(disp))/
		disp+n*lgamma(1/disp)+sum(log(y)*wt),
	"inverse Gauss"=maxlike <- (sum(wt)*(log(disp*2*pi)+1)+3*sum(log(y)*wt))/2)
	maxlike <- maxlike-sum(log(delta))
	fitted.values <-  as.vector(mu(z0$estimate))
	residuals <-  y-fitted.values
	if(np==1)cov <- 1/z0$hessian
	else {
		a <- qr(z0$hessian)
		if(a$rank==np)cov <- solve(z0$hessian)
		else cov <- matrix(NA,ncol=np,nrow=np)}
	cov <- 2*cov*z0$minimum/sum(wt)
	se <- sqrt(diag(cov))
	z1 <- list(
		call=call,
		dist=dist,
		delta=delta,
		mu=mu,
		prior.weights=wt,
		maxlike=maxlike,
		dispersion=disp,
		fitted.values=fitted.values,
		residuals=residuals,
		aic=maxlike+np+1,
		df=sum(wt)-np,
		coefficients=z0$estimate,
		np=np,
		se=se,
		cov=cov,
		corr=cov/(se%o%se),
		gradient=z0$gradient,
		iterations=z0$iterations,
		code=z0$code)
	class(z1) <- "nlr"
	return(z1)}

residuals.nlr <- function(z) z$residuals
fitted.values.nlr <- function(z) z$fitted.values
coefficients.nlr <- function(z) z$coefficients
weights.nlr <- function(z) z$prior.weights
df.residual.nlr <- function(z) z$df
deviance.nlr <- function(z) 2*z$maxlike

print.nlr <- function(z) {
	cat("\nCall:\n",deparse(z$call),"\n\n",sep="")
	if(z$code>2)cat("Warning: no convergence - error",z$code,"\n\n")
	cat(z$dist,"distribution\n\n")
	t <- deparse(z$mu)
	cat("Mean function:",t[2:length(t)],sep="\n")
	cat("\n-Log likelihood   ",z$maxlike,"\n")
	cat("Degrees of freedom",z$df,"\n")
	cat("AIC               ",z$aic,"\n")
	cat("Iterations        ",z$iterations,"\n\n")
	cat("Mean parameters:\n")
	coef.table <- cbind(z$coefficients[1:z$np], z$se[1:z$np])
	dimnames(coef.table) <- list(seq(1,z$np), c("estimate", "se"))
	print.default(coef.table, digits=4, print.gap=2)
	cat("\nDispersion estimate:",z$dispersion,"\n")
	if(z$np>1){
		cat("\nCorrelations:\n")
		dimnames(z$corr) <- list(seq(1,z$np),seq(1,z$np))
		print.default(z$corr, digits=4)}
	invisible(z)}
nordr <- function(y, dist="proportional", mu, linear=NULL, pmu, 
      pblock, wt=NULL, print.level=0, ndigit=10, gradtol=0.00001,
      steptol=0.00001, fscale=1, iterlim=100, typsiz=abs(p),
      stepmax=10*sqrt(p%*%p)){
lf <- function(p){
   g <- exp(mu1(p[1:npl])+block%*%p[npl1:np])
   g <- g/(1+g)
   if(mdl==1){
	g <- c(g,ext)
	g <- g[1:nlen]/g[nrows1:nlenr]}
   -sum(pwt*(resp*log(g)+(1-resp)*log(1-g)))
}
lf3 <- function(p){
    mu <- mu1(p[1:npl])
    g <- exp(mu*(y-1)+resp%*%p[npl1:np])/
      exp(mu%o%0:my+matrix(rep(cumsum(c(0,0,p[npl1:np])),nrows),ncol=my+1,byrow=T))%*%ext
    -sum(pwt*log(g))}
call <- sys.call()
dist <- match.arg(dist,c("proportional odds","continuation ratio",
     "adjacent categories"))
if(min(y)!=1)stop("ordinal values must be start at 1")
else if(any(y!=trunc(y)))stop("ordinal values must be integers")
else my <- max(y)-1
if(dist=="proportional odds")mdl <- 1
else if(dist=="continuation ratio")mdl <- 2
else mdl <- 3
nrows <- length(y)
nrows1 <- nrows+1
nlen <- my*nrows
nlenr <- nlen+nrows
nlp <- npl <- length(pmu)
npl1 <- npl+1
if(is.language(mu))linear <- mu
if(is.language(linear)){
	mt <- terms(linear)
	if(is.numeric(mt[[2]])){
		if(!is.function(mu)){
			mu1 <- function(p) p[1]*rep(1,nrows)
			nlp <- 1}
		else mu1 <- function(p) mu(p, p[1]*rep(1,nrows))}
	else {
		mf <- model.frame(mt,sys.frame(sys.parent()),na.action=na.fail)
		dm1 <- model.matrix(mt, mf)
		npt1 <- ncol(dm1)
		if(!is.function(mu)){
			mu1 <- function(p) as.vector(dm1%*%p[1:npt1])
			nlp <- npt1}
		else mu1 <- function(p) mu(p, dm1%*%p[1:npt1])}}
else if(!is.function(mu)){
	mu1 <- function(p) p[1]*rep(1,nrows)
	nlp <- 1}
else {
	mu1 <- mu
	if(length(mu1(pmu))==1)mu1 <- function(p) mu(p)*rep(1,nrows)}
if(nlp!=npl)
	stop("Number of initial estimates for mu does not correspond to model")
if(mdl==1)ext <- rep(1,nrows)
else if(mdl==3)ext <- rep(1,my+1)
if(mdl==3)resp <- NULL
else resp <- matrix(as.integer(y==1),ncol=1)
block <- NULL
pwt <- matrix(as.integer(y<3),ncol=1,nrow=nrows)
for(i in 2:my){
      resp <- cbind(resp,as.integer(y<=i))
      block <- cbind(block,as.integer(c(rep(0,nrows*(i-1)),
	    rep(1,nrows),rep(0,nrows*(my-i)))))
      pwt <- cbind(pwt,as.integer(y<i+2))}
if(mdl!=1)resp <- 1-resp
if(mdl!=3)resp <- as.vector(resp)
pwt <- as.vector(pwt)
if(!is.null(wt)){
	if(!is.vector(wt))stop("wt must be a vector")
	else if(length(wt)!=nrows)stop("wt must have length",nrows)
	if(mdl==3)pwt <- wt
	else pwt <- rep(wt,my)*pwt}
if(missing(pblock)||length(pblock)!=my-1)
	stop(paste(my-1,"initial values of block parameters must be supplied"))
p <- c(pmu,pblock)
np <- length(p)
if(mdl==3)z <- nlm(lf3, p, hessian=T, print.level=print.level,
	typsiz=typsiz, ndigit=ndigit, gradtol=gradtol, stepmax=stepmax,
	steptol=steptol, iterlim=iterlim, fscale=fscale)
else z <- nlm(lf, p, hessian=T, print.level=print.level,
	typsiz=typsiz, ndigit=ndigit, gradtol=gradtol, stepmax=stepmax,
	steptol=steptol, iterlim=iterlim, fscale=fscale)
maxlike <- z$minimum
if(!is.language(linear))cname <- paste("p",1:npl,sep="")
else {
     cname <- colnames(dm1)
     if(is.function(mu))
	cname <- c(cname,paste("p",(length(cname)+1):npl,sep=""))}
a <- qr(z$hessian)
if(a$rank==np)cov <- solve(z$hessian)
else cov <- matrix(NA,ncol=np,nrow=np)
se <- sqrt(diag(cov))
corr <- cov/(se%o%se)
dimnames(corr) <- list(1:np,1:np)
z1 <- list(
   call=call,
   dist=dist,
   wt=wt,
   maxlike=maxlike,
   aic=maxlike+np,
   mu=mu,
   linear=linear,
   coefficients=z$estimate[1:npl],
   cname=cname,
   np=np,
   npl=npl1-1,
   nrows=nrows,
   block=z$estimate[npl1:np],
   cov=cov,
   corr=corr,
   se=se,
   iterations=z$iter,
   code=z$code)
class(z1) <- "nordr"
z1}

print.nordr <- function(z, digits = max(3, .Options$digits - 3)){
	m <- z$states
	cat("\nCall:\n",deparse(z$call),"\n\n",sep="")
	cat(z$dist,"model\n\n")
	if(z$code>2)cat("Warning: no convergence - error",z$code,"\n\n")
	cat("-Log likelihood   ",z$maxlike,"\n")
	cat("AIC               ",z$aic,"\n")
	cat("Iterations        ",z$iterations,"\n")
	cat("\nMean coefficients\n")
	if(is.function(z$mu)){
		t <- deparse(z$mu)
		cat("Mean function:\n",t[2:length(t)],sep="\n")
		if(is.language(z$linear))
			cat("Linear part: ",deparse(z$linear[[2]]),"\n")}
	coef.table <- cbind(z$coef,z$se[1:z$npl])
	dimnames(coef.table) <- list(z$cname,c("estimate","s.e."))
	print.default(coef.table, digits=digits, print.gap=2)
	cat("\nBlock coefficients\n")
	coef.table <- cbind(z$block,z$se[(z$npl+1):z$np])
	dimnames(coef.table) <- list(paste("b",2:(z$np-z$npl+1),sep=""),
			     c("estimate","s.e."))
	print.default(coef.table, digits=digits, print.gap=2)
	cat("\nCorrelation matrix\n")
	print.default(z$corr, digits=digits)
}
# standard pharmacokinetic models
#
# open zero-order one-compartment model
# p[1]: log volume (V)
# p[2]: log elimination rate (ke)
# end:  time when infusion stops
mu1.0o1c <- function(p, times, dose=1, end=0.5) {
	ke <- exp(p[2])
	dose/(exp(p[1])*ke)*((1-exp(-ke*times))*(times<=end)+
		(1-exp(-ke*end))*exp(-ke*(times-end))*(times>end))}
#
# open first-order one-compartment model
# p[1]: log volume (V)
# p[2]: log absorption rate (ka)
# p[3]: log elimination rate (ke)
mu1.1o1c <- function(p, times, dose=1) {
	ka <- exp(p[2])
	ke <- exp(p[3])
	exp(p[2]-p[1])*dose/(ka-ke)*(exp(-ke*times)-exp(-ka*times))}
#
# open first-order two-compartment model (ordered)
# p[1]: log volume (V)
# p[2]: log absorption rate (ka)
# p[3]: log elimination rate (ke)
# p[4]: log transfer rate between compartments (k12)
mu1.1o2c <- function(p, times, dose=1) {
	ka <- exp(p[2])
	ke <- exp(p[3])
	k12 <- exp(p[4])
	ka*k12*exp(-p[1])*dose/(k12-ka)*((exp(-ka*times)-exp(-ke*times))/
		(ke-ka)-(exp(-k12*times)-exp(-ke*times))/(ke-k12))}
#
# open first-order two-compartment model (ordered, absorption and transfer equal)
# p[1]: log volume (V)
# p[2]: log absorption rate (ka)
# p[3]: log elimination rate (ke)
mu1.1o2cl <- function(p, times, dose=1) {
	ka <- exp(p[2])
	ke <- exp(p[3])
	ka^2*exp(-p[1])*dose/(ka-ke)*((exp(-ka*times)-exp(-ke*times))/(ke-ka)
		-times*exp(-ka*times))}
#
# open first-order two-compartment model (circular)
# p[1]: log volume (V)
# p[2]: log absorption rate (ka)
# p[3]: log elimination rate (ke)
# p[4]: log rate to second compartment (k12)
# p[5]: log rate from second compartment (k21)
mu1.1o2cc <- function(p, times, dose=1) {
	ka <- exp(p[2])
	ke <- exp(p[3])
	k12 <- exp(p[4])
	k21 <- exp(p[5])
	beta <- (k12+k21+ke-sqrt((k12+k21+ke)^2-4*k21*ke))/2
	alpha <- (k21*ke)/beta
	exp(p[2]-p[1])*dose*((k21-alpha)*exp(-alpha*times)/
		((ka-alpha)*(beta-alpha))+(k21-beta)*exp(-beta*times)/
		((ka-beta)*(alpha-beta))+(k21-ka)*exp(-ka*times)/
		((beta-ka)*(beta-ka)))}
#
# simultaneous models for parent drug and metabolite
#
# zero-order one-compartment model
# p[1]: log parent drug volume (Vp)
# p[2]: log parent drug direct elimination rate (kep)
# p[3]: log transformation rate from parent to metabolite (kpm)
# p[4]: log metabolite elimination rate (kem)
# p[5]: log metabolite volume (Vm)
# ind:  indicator vector - 1 for parent, 0 for metabolite
# end:  time when infusion stops
mu2.0o1c <- function(p, times, dose=1, ind, end=0.5) {
	Vp <- exp(p[1])
	kpm <- exp(p[3])
	kp <- exp(p[2])+kpm
	kem <- exp(p[4])
	Vm <- exp(p[5])
	kemp <- kem-kp
	tmp1 <- exp(-kp*times)
	tmp2 <- kpm/(kp*kem*Vm)
	g1 <- exp(-kp*end)
	g2 <- exp(-kem*end)
	cend <- (1-g1)/(Vp*kp)
	cexp <- exp(-kp*(times-end))*(times>end)
	cmend <- tmp2*(1+kp/kemp*g2-kem/kemp*g1)
	tmp3 <- cend*kpm*Vp/(kemp*Vm)
	dose*(ind*((1-tmp1)/(Vp*kp)*(times<=end)+cend*cexp)+
		(1-ind)*(tmp2*(1+kp/kemp*exp(-kem*times)-kem/kemp*
		tmp1)*(times<=end)+(g2/g1*cexp*tmp3+
		(cmend-tmp3)*exp(-kem*(times-end))/g2*(times>end))))}
#
# zero-order two-compartment model
# p[1]: log parent drug volume (Vp)
# p[2]: log parent drug direct elimination rate (kep)
# p[3]: log parent drug rate to second compartment (k12)
# p[4]: log parent drug rate from second compartment (k21)
# p[5]: log transformation rate from parent to metabolite (kpm)
# p[6]: log metabolite elimination rate (kem)
# ind:  indicator vector - 1 for parent, 0 for metabolite
# end:  time when infusion stops
mu2.0o2c <- function(p, times, dose=1, ind, end=0.5) {
	Vp <- exp(p[1])
	k12 <- exp(p[3])
	k21 <- exp(p[4])
	kpm <- exp(p[5])
	kp <- exp(p[2])+kpm
	kem <- exp(p[6])
	tmp1 <- sqrt((kp+k12+k21)^2-4*k21*kp)
	lam1 <- 0.5*(kp+k12+k21+tmp1)
	lam2 <- 0.5*(kp+k12+k21-tmp1)
	tmp2 <- (1-exp(-kem*times))/kem
	tmp3 <- (1-exp(-kem*end))/kem
	tmp4 <- lam1-k21
	tmp5 <- lam2-k21
	tmp6 <- exp(-lam1*times)
	tmp7 <- exp(-lam2*times)
	tmp8 <- exp(-lam1*end)
	tmp9 <- exp(-lam2*end)
	tmp10 <- exp(-kem*times)
	tmp11 <- tmp6/tmp8
	tmp12 <- tmp7/tmp9
	tmp13 <- exp(-kem*end)
	tmp14 <- exp(-kem*(times-end))
	dose/Vp*(ind*((tmp4*(1-tmp6)/(lam1*tmp1)-tmp5*(1-tmp7)/(lam2*tmp1))*
		(times<=end)+(tmp4*(1-tmp8)*tmp11/(lam1*tmp1)-tmp5*(1-tmp9)*
		tmp12/(lam2*tmp1))*(times>end))+(1-ind)*kpm*(((tmp4*tmp2-(tmp6-
		tmp10)/(kem-lam1))/(lam1*tmp1)-(tmp5*tmp2-(tmp7-tmp10)/
		(kem-lam2))/(lam1*tmp1))*(times<=end)+(((tmp4*tmp3-(tmp8-
		tmp13)/(kem-lam1))/(lam1*tmp1)-(tmp5*tmp3-(tmp9-tmp13)/
		(kem-lam2))/(lam1*tmp1))*tmp14+tmp4*(1-tmp8)*
		(tmp14-tmp11)/(lam1*(lam1-kem)*tmp1)-tmp5*(1-tmp9)*
		(tmp14-tmp12)/(lam2*(lam2-kem)*tmp1))*(times>end)))}
#
# first-order one-compartment model
# p[1]: log volume (V)
# p[2]: log parent drug absorption rate (kap)
# p[3]: log parent drug direct elimination rate (kep)
# p[4]: log transformation rate from parent to metabolite (kpm)
# p[5]: log metabolite elimination rate (kem)
# ind:  indicator vector - 1 for parent, 0 for metabolite
mu2.1o1c <- function(p, times, dose=1, ind) {
	kap <- exp(p[2])
	kep <- exp(p[3])
	kem <- exp(p[5])
	kap*exp(p[1])*dose/(kap-kep)*(ind*(exp(-kep*times)-exp(-kap*times))+
	(1-ind)*exp(p[4])*(exp(-kap*times)/(kap-kem)-exp(-kep*times)/(kep-kem)+
	(1/(kep-kem)-1/(kap-kem))*exp(-kem*times)))}
#
# zero-order one-compartment first-pass model
# p[1]: log parent drug volume (Vp)
# p[2]: log parent drug direct elimination rate (kep)
# p[3]: log transformation rate from parent to metabolite (kpm)
# p[4]: log metabolite elimination rate (kem)
# p[5]: log metabolite volume (Vm)
# p[7]: logit of proportion going to first pass (pfp)
# ind:  indicator vector - 1 for parent, 0 for metabolite
# end:  time when infusion stops
mu2.0o1cfp <- function(p, times, dose=1, ind, end=0.5) {
	Vp <- exp(p[1])
	kpm <- exp(p[3])
	kp <- exp(p[2])+kpm
	kem <- exp(p[4])
	Vm <- exp(p[5])
	kemp <- kem-kp
	tmp1 <- exp(-kp*times)
	tmp2 <- exp(-kem*times)
	tmp3 <- kpm/(kp*kem*Vm)
	g1 <- exp(-kp*end)
	g2 <- exp(-kem*end)
	cend <- (1-g1)/(Vp*kp)
	cexp <- exp(-kp*(times-end))*(times>end)
	cmend <- tmp3*(1+kp/kemp*g2-kem/kemp*g1)
	tmp4 <- cend*kpm*Vp/(kemp*Vm)
	pfp <- exp(p[6])
	lpfp <- pfp/(1+pfp)
	dose*(ind*((1-tmp1)/(Vp*kp)*(times<=end)+cend*cexp)*lpfp+
		(1-ind)*(((1-tmp2)/(Vm*kem)*(times<=0.5)+(1-g2)/(Vm*kem)*
		exp(-kem*(times-0.5))*(times>0.5))/(1+pfp)+(tmp3*
		(1+kp/kemp*tmp2-kem/kemp*tmp1)*(times<=end)+(g2/g1*cexp*tmp4+
		(cmend-tmp4)*exp(-kem*(times-end))/g2*(times>end)))*lpfp))}
#
# first-order one-compartment first-pass model
# p[1]: log volume (V)
# p[2]: log parent drug absorption rate (kap)
# p[3]: log parent drug direct elimination rate (kep)
# p[4]: log transformation rate from parent to metabolite (kpm)
# p[5]: log metabolite first-pass elimination rate (kefp)
# p[6]: log metabolite elimination rate (kem)
# p[7]: logit of proportion going to first pass (pfp)
# ind:  indicator vector - 1 for parent, 0 for metabolite
mu2.1o1cfp <- function(p, times, dose=1, ind) {
	kap <- exp(p[2])
	kep <- exp(p[3])
	kefp <- exp(p[5])
	kem <- exp(p[6])
	pfp <- exp(p[7])
	lpfp <- pfp/(1+pfp)
	exp(p[2]-p[1])*dose*(ind*(exp(-kep*times)-exp(-kap*times))/
		(kap-kep)*lpfp+(1-ind)*((exp(-kefp*times)-exp(-kap*times))/
		(kap-kefp)/(1+pfp)+exp(p[4])*(exp(-kap*times)/(kap-kem)-
		exp(-kep*times)/(kep-kem)+(1/(kep-kem)-1/(kap-kem))*
		exp(-kem*times))/(kap-kep)*lpfp))}
rs2 <- function(y,x1,x2,power=c(1,1),weight=rep(1,length(x1)),family=normal){
	if(length(power)!=2){
		cat("Two estimates of power parameters must be supplied\n")
		return(1)}
	a <- power[1]
	b <- power[2]
	test <- T
	while(test){
		xx1 <- x1^a
		xx2 <- x2^b
		u <- glm(y~xx1+xx2+I(xx1^2)+I(xx2^2)+xx1:xx2,family=family,
			weight=weight)
		z1 <- (u$coef[2]*xx1+2*u$coef[4]*xx1^2+u$coef[6]*xx1*xx2)*
			log(x1)
		z2 <- (u$coef[3]*xx2+2*u$coef[5]*xx2^2+u$coef[6]*xx1*xx2)*
			log(x2)
		u <- glm(y~xx1+xx2+I(xx1^2)+I(xx2^2)+xx1:xx2+z1+z2,
			family=family,weight=weight)
		a <- a+u$coef[6]
		b <- b+u$coef[7]
		test <- (u$coef[6]^2>0.00001)||(u$coef[7]^2>0.00001)}
	u <- glm(y~xx1+xx2+I(xx1^2)+I(xx2^2)+xx1:xx2,family=family,
		weight=weight)
	u$df.residual <- u$df.residual-2
	z <- list(model=u,estimates=c(a,b))
	class(z) <- "rs"
	return(z)}

rs3 <- function(y,x1,x2,x3,power=c(1,1,1),weight=rep(1,length(x1)),
	family=normal){
	if(length(power)!=3){
		cat("Three estimates of power parameters must be supplied\n")
		return(1)}
	a <- power[1]
	b <- power[2]
	d <- power[3]
	test <- T
	while(test){
		xx1 <- x1^a
		xx2 <- x2^b
		xx3 <- x3^d
		xx12 <- xx1*xx2
		xx13 <- xx1*xx3
		xx23 <- xx2*xx3
		u <- glm(y~xx1+xx2+xx3+I(xx1^2)+I(xx2^2)+I(xx3^2)+
			xx12+xx13+xx23,family=family,weight=weight)
		z1 <- (u$coef[2]*xx1+2*u$coef[5]*xx1^2+u$coef[8]*xx12+
			+u$coef[9]*xx13)*log(x1)
		z2 <- (u$coef[3]*xx2+2*u$coef[6]*xx2^2+u$coef[8]*xx12+
			u$coef[10]*xx23)*log(x2)
		z3 <- (u$coef[4]*xx2+2*u$coef[7]*xx2^2+u$coef[9]*xx13+
			u$coef[10]*xx23)*log(x3)
		u <- glm(y~xx1+xx2+xx3+I(xx1^2)+I(xx2^2)+I(xx3^2)+
			xx12+xx13+xx23+z1+z2+z3,family=family,weight=weight)
		a <- a+u$coef[11]
		b <- b+u$coef[12]
		d <- d+u$coef[13]
		test <- (u$coef[11]^2>0.00001)||(u$coef[12]^2>0.00001)||
			(u$coef[13]^2>0.00001)}
	u <- glm(y~xx1+xx2+xx3+I(xx1^2)+I(xx2^2)+I(xx3^2)+xx12+xx13+xx23,
		family=family,weight=weight)
	u$df.residual <- u$df.residual-3
	z <- list(model=u,estimates=c(a,b,d))
	class(z) <- "rs"
	invisible(z)}

print.rs <- function(z){
	cat("Powered transformed response surface\n\n")
	cat("Powers:",z$est,"\n")
	cat("AIC   :",z$model$dev+2*(length(z$model$coef)+length(z$est)),"\n\n")
	print(summary(z$model,corr=F))}

