#            util.R Utilty routines for ORNL/TM-2004/146 
#
#     Statistical Analyis of Data with Non-Detectable Values
#              E. L. Frome and J. P. Watkins
#     Revision 6: 28 July 2004
# 
#  Name       Purpose
#  plend()    Product Limit Estimate of CDF for data with non-detects
#  plquan()   calculate the qth quantile from PLE of CDF
#  lnclxpnd() calculate confidence intervals for lognormal percentiles 
#  kmms()     calculate Kaplan-Meier mean and CLs
#  coxcl()    calculate modified Cox Interval for lognormal mean
#  LKcl()     calculate CLs for lognormal mean: Lyles and Kupper 
#  lnstats()  calculate summary statistic for left censored sample 
#  qqlognA()  q-q plot for left censored lognormal sample
#  lnexh1     Exhibit 1 from ORNL/TM-2004/146
#  lnexh2     Exhibit 2 from ORNL/TM-2004/146
#  readss()   read data in spread sheet format
#
############################ plend ####################################

plend<-function(dd=ex1){
#    Product Limit Estimate  for positive data with  non-detectable
#   ( left censored data) values see Section 3.5
# USAGE: plend(dd)
# ARGUMENTS: dd is an n by 2 matrix or data frame
#       dd[,1]=  exposure variable in column 1
#       dd[,2] = censor ( 0 for non-detect 1 for detect)
# VALUE:   data frame with columns
#          apl    a   ple  n r   surv
#
#      apl[j] is adjusted ple used in q-q plot
#      a(j)   is the  value of jth detect (ordered)
#      ple(j) is product limit estimate of cdf (Kaplin-Meier)
#      n(j) = number of detects or non-detects <= x(j)
#      r(j)=  number of detects at x(j)
#      suv(j) = 1 - ple(j)
# REFERENCES: Schmoyer et al (1996) Environmental and Ecological
#             Statistics,3 81-79
# NOTE: see R package dblcens
nn<- length(dd[,1]); n<-rev( 1: nn )
t1<- as.matrix( cbind( dd[ rev(order(dd[,1]) ),1:2],n ) )
dimnames(t1)<- list(NULL,c("x","cen","n") )
t2<- t1[ t1[,2]==1,]  # select detects
ung<- c(-1, diff(t2[,1] )) 
d<- rev( table( t2[,1] ) ) 
#  columns of t2 are x[j] n[j] d[j] for cen=1
#   IN REVERSE order
t2<- as.matrix( cbind( t2[ ung< 0,],d )[,c(1,3:4)] )
#  x contains detects and nx is number of detects
x<- dd[ dd[,2]==1,1] ; nx<-length(x) 
ndt<- max(x)
# add row if min(x) is a non-detect
if( nx < dim(dd)[1] ){
ndt<- dd[ dd[,2]==0, 1] #  nondetects
if( min(ndt) < min(x) ){
nl<- sum( ifelse(dd[,1] <= min(ndt),1,0) )
 frow<- c(min(ndt),nl,0)
 t2<- rbind(t2,frow)
}}
# ple = Prod( [n[j] - d[j] ]/n[j] ) 
n1<- dim(t2)[1]
ple0<- cumprod( (t2[,2]-t2[,3])/t2[,2] )
ple<- c(1,ple0[1:n1-1])
apl<- as.matrix( rev( (ple0 + ple)/2 ) )
#dimnames(apl)<-list( NULL,c("apl","a","ple","n", "r") )
suv<- 1 - ple ;  suv<-c(suv[2:n1],1 )
t2<- cbind(t2,ple,suv)
t2<- t2[order(t2[,1]),]
#t2<- data.frame(apl,t2[,c(1,4,2,3,5)])
t2<- cbind(apl,t2[,c(1,4,2,3,5)])
dimnames(t2)<-list( NULL,c("apl","a","ple","n", "r","surv") )
t2<-data.frame(t2)
t2
}
############################ plquan ####################################

plquan<-function(pe=ple,qq=0.95){
# Find the qth quantile from PLE of CDF
# USAGE: plend(pe=ple,qq=0.95)
# ARGUMENTS: pe is data.frame from plend
#            qq is such that qq*100 is percentile
# VALUE: Xq the 100qqth percentile 
x <- c(0,pe[,2])
p <- c(0,pe$ple)
j <- 0
for( jj in 1:length(x) ) {j<- j+1
   if( p[jj] > qq)  break}
 xq<- (x[j] - x[j-1])/(p[j] -p[j-1])
 xq<- x[j-1] + xq*( qq - p[j-1])
xq
}

############################ lnclxpnd ####################################

lnclxpnd<-function(ae=met,p=95,gam=95){
#  calculate with confidence level gam that 
#  p percent of population lies  below the upper bound xpu
#     and above the lower bound xpl so that
#  (xpl,xpu) is an approximate 100*[ 1 - 2*(100-p)/100 ] percent  
#  Confidence Interval for the Xp percentile of lognormal distribution
# USAGE: lnclxpnd(ae,p,gam)
# ARGUMENTs:ax is output from mlndln()
#           p is percentile and 
#           gam is confidence level
# VALUE: Xp is pth percentile of lognormal distribution
#        (xpl,xpu) METHOD 1 Confidence Interval for Xp
#        (expl,expu) METHOD 2 Confidence Interval for Xp
# NOTE: The upper bound xpu is the UTL ( upper tolerance limit)
#   for the pth percentile, i.e UTL-pg
#  METHOD 1 IS LARGE SAMPLE RESULT
#  x is lognormal so y=log(x) is normal(mean=ym,sd=ysd)
#    yp = ym + zp*ysd is pth percentile of y
#    ypu = ym + t(gam/100,nt1-1 )*zp*sdyq  is upper bound
#  var(yqu) = var(yb) + zp^2*var(ysd) +2*zp*cov(yb,ysd)
#   zp is pth percentile of N(0,1) and t() is quantile of  t distn
# METHOD 2 uses tolerance limit factor fron non-central t distribution
#  these are exact limits when there are no nondetects
# REQUIRESs extol()
mu <- ae[1,1] ; sig <- ae[1,2] ; m <-ae[2,6] 
vmu <- ae[2,1]^2 ;  vsig <- ae[2,2]^2
zt <- qnorm(p/100) ; yp<- mu + zt*sig
#  cov(mu,sig) from mlndln() is in ae[2,5]
sdyp <- sqrt( vmu + zt^2*vsig + 2*zt*ae[2,5] )
xp <- exp(yp)
xpu <- exp( yp + qt(gam/100,m-1 )*sdyp )
xpl <- exp( yp - qt(gam/100,m-1 )*sdyp )
#  km <- (log(xpu) -mu)/sig  # km is equivalent to K factor
K<-extol(m,p/100,gam/100)[1] # apprx K value
expu <- exp(mu + K*sig)  
expl <- exp(mu + sig*extol(m,p/100,(100 -gam)/100)[1] ) 
#out<-c(xp,xpl,xpu,expl,expu,km,K,m)
#names(out)<-c("Xp","xpl","xpu","expl","expu","km","K","m")
out<-c(xp,xpl,xpu,expl,expu)
names(out)<-c("Xp","xpl","xpu","expl","expu")
out
}

############################ kmms  ####################################

kmms <- function(dd=ex1,cl=95){
#    Kaplan- Meier(KM) estimate of mean and Standard Error of
#    the mean for left censored data  see Section 3.5
# USAGE: kmms(dd,cl)
# ARGUMENTS:  dd[,1]= data(exposure variable)
#          dd[,2] = censor ( 0 for non-detect 1 for detect)
#          cl= one-sided confidence level percent
# VALUE:
#  KM-mean = KM non-parametric estimate of mean 
#  KM-se = standard error of mean (adjusted)
#  KM-L and KM-U lower and upper confidence limits 
# NOTE: (kml,kmu) is a 100*[1 - (1-cl/100)*2] Percent Confidence Interval
# REFERENCES:
#  Kaplan- Meier (1958) J Am Stat Assoc 457-481
#       Turnbull (1976) J Royal Stat Soc 290-295
nn<- length(dd[,1]); n<-rev( 1: nn )
#  t1 columns of x cen n in reverse order
t1<- cbind( dd[ rev(order(dd[,1]) ),1:2],n )
t2<- t1[ t1[,2]==1,]  # select detects
ung<- c(-1, diff(t2[,1] )) 
d<- as.vector( rev( table( t2[,1] ) ) )
# t2 keep rows with d GT O 
t2<- as.matrix( cbind( t2[ ung< 0,],d )[,c(1,3:4)] )
#  x contains detects and nx is number of detects
x<- dd[ dd[,2]==1,1] ; nx<-length(x) 
ndt<-max(x)
#  if min val of x is non-detect add a row
if( nx < dim(dd)[1] ){
ndt<- dd[ dd[,2]==0, 1] #  nondetects
if( min(ndt) < min(x) ){
nl<- sum( ifelse(dd[,1] <= min(ndt),1,0) )
 frow<-c(min(ndt),nl,0)
 t2<-rbind(t2,frow)
}}
# caculate product limit estimate
  ple<-cumprod( (t2[,2]-t2[,3])/t2[,2] )
if( min(ndt) < min(x) ) ple[ dim(t2)[1] ] <- 0 
sv<- 1 - ple
xv<- abs(diff(c(t2[,1],0)) )
a<- sv*xv  # a is area for jth interval
ac<- xv*ple ;n1<-dim(t2)[1]
  # if( t2$d[nl==0 )ac[n1-1]<-0
B<- rev(cumsum(rev(ac)))
t2<- cbind(t2,ple,sv,a,ac,B)
t2<- data.frame( t2[order(t2[,1]),] )
#  compute terms for variance
vb<-ifelse( (t2[,2] -t2[,3])==0 , 0, 1/((t2[,2] -t2[,3])*t2[,2]) )
vb<- t2$d*t2$B^2*vb
t2<-data.frame(t2,vb)
kmm<-sum( t2$a) # Kaplan-Meier mean
kmvb<-sum(vb)   # Kaplan-Meier variance (unadjusted)
kmseu<- sqrt(kmvb); kmse<- kmseu* sqrt(nx/(nx-1) )
# calculate two-sided cl percent CLs for mean
# use qt( 1 - ( 1 - cl/100)/2  ,nx-1
# 
cd<- qt( cl/100 ,nx-1)*kmse
kml<- kmm - cd ; kmu<- kmm + cd
stat<-c(kmm,kmse,kml,kmu,cl)
names(stat)<-c("KM-mean","KM-se","KM-L","KM-U","CL-one")
stat
}

############################ coxcl  ####################################

coxcl<-function(du=ex1,cl=95,dat=T){
#    Calculate confidence limits for lognormal mean
#    using equivalent of Cox's direct method modified 
#    for non-detects as described in Section 3.2
# USAGE: coxcl(du,cl,dat)
# ARGUMENTS: if dat=T  du is an n by 2 matrix or data frame
#            if dat=F  du is matrix from mlndln(du)
#            cl is one-sided confidence level(%) 
# VALUE: ML estimate of AM= exp(mu +0.5sig^2)
#        AM-L lower confidence limit for AM
#        AM-U upper confidence limit for AM
# NOTE: (AM-L,AM-U) is a 100*[1 - (1-cl/100)*2] Percent Confidence Interval
    if( dat==F ) ys <- du
    if( dat==T ) ys<- mlndln(du)
#   see mlndln() comments ys[1,3] is ML estimate of logED
#    ys[2,3] is standard of logED
m <- ys[2,6]  # m = number of non-detects
lex <- ys[1,3]  # same as mu + 0.5*sig^2
tv <- qt(cl/100, m-1)
lexl <- lex - tv*ys[2,3]
lexu <- lex + tv*ys[2,3]
out<- c( exp( c(lex,lexl,lexu) ) )
names(out)<-c("AM","AM-L","AM-U")
out
}
############################ LKcl ####################################

LKcl<-function(din=ex1,cl=95,dat=T){
#  Calculate confidence limits for lognormal mean
#  using approximate method described in Section 3.2
# USAGE: LKcl(din,cl,dat)
# ARGUMENTS: if dat=T  din is an n by 2 matrix or data frame
#            if dat=F  du is matrix from mlndln(din)
#            cl is one-sided confidence level(%)
# VALUE: ML estimate of AM= exp(mu +0.5sig^2)
#        LK-L lower confidence limit for AM
#        LK-U upper confidence limit for AM
# NOTE: (LK-L,LK-U) is a 100*[1 - (1-cl/100)*2] Percent Confidence Interval
# If there are no non-detects (m=n) then (LK-L,LK-U) will be
# an approximation to Land's exact method
# REFERENCES:
#     Lyles and Kupper (1996) JAIHA vol 57 6-15
ml <- din
if( dat==T ) ml <- mlndln(din) 
  yb<-ml[1,1] ; sy<-ml[1,2]; vyb<-ml[2,1]^2
#  
  n <- ml[2,6]  # number of non-detects
 alp <- 1 -cl/100 
 ch<- sqrt( qchisq(alp,n-1) ) ; clo<-sqrt( qchisq(1-alp,n-1) )
#  see Section 3.2 equation 11
cuh<- (sqrt(n)*(sy/2)*sqrt((n-1)/n))/ch  + qt(1-alp,n-1)/sqrt(n)
cul<- (sqrt(n)*(sy/2)*sqrt((n-1)/n))/clo  + qt(alp,n-1)/sqrt(n)
lex<- yb +0.5*sy^2
lexl<-yb + cul*sy ; kl<-  sqrt(n-1)*(lexl - lex)/sy
lexu<- yb + cuh*sy ; ku<- sqrt(n-1)*(lexu - lex)/sy
out<-c(exp(lex) ,exp(lexl),exp(lexu) )
names(out)<-c("AM","LK-L","LK-U")
out
}

############################ lnstats  ####################################
lnstats <- function(d=ex1,OEL=1250,pc=95,cl=95,mth=1){
# lnstats.R 17 SEP 2003 rev 10 Feb 2004 rev 7 Jul 04
#  calculate summary statistic for left censored sample x[]
#  from lognormal distribution based on Maximum Likelihood
#  and nonparametric methods
#  INPUT:
#    d is 2col data frame or matrix  with
#    x  in column 1  and cen=1 for  detect 0 for nondetect column 2
#    cl is confidence level (one sided) for confidence intervals
#    pc is quantile for UTL-pc-cl
#    mth method =1 for asymptotic =2 ("exact" method)
# REQUIRES mlndln(d) to calculate Maximum Likelihood Estimates
#   plend()  product limit estimate of CDF for left censored data
#

nt <- length(d[,1])
ndet <- sum(d[,2])
PCndet <- round( 100*(nt-ndet)/nt,1)
du <- d[,1:2] 

ys <- mlndln(du) # row 1 and 2 bias adjusted
GM <- exp(ys[1,1])
GSD <- exp(ys[1,2])

#  Method 1 is equivalent to modified  Cox direct Method
if(mth==1) {clo<-coxcl(ys,cl,F)
   EX<-clo[1]
   EXL<-clo[2]; nEXL=paste("LCLc",pc,sep="-")
   EXU<- clo[3]; nEXU=paste("UCLc",pc,sep="-")
}

#  If mth=2 Use Modifed Method of Lyles and Kupper 
#  to calculate approximate CLs for lognormal mean
if(mth==2) {lko<-LKcl(ys,cl,F)
   EX<-lko[1]
   EXL<-lko[2]; nEXL=paste("LCLa",pc,sep="-")
   EXU<- lko[3]; nEXU=paste("UCLa",pc,sep="-")
}
#  calculate Product limit estimate of CDF for left censored data
cdf <- plend(du)
#  RH is (progressively left censored version of correlation coef)
#  see Verrill and Johnson JASA 1988 Equation 4.4
#  in no censoring Rsq= RH^2 will be approx equal to Shapiro-Wilk W
Rsq <- cor(qnorm(cdf[,1]),log(cdf[,2]) )^2

xtmp<-lnclxpnd(ys,p=pc,gam=cl)
pc.obs<- plquan( cdf  , pc/100)	# find pcth percentile of x
pc.est<-xtmp[1]

pco<-paste("Obs%",pc,sep=""); pce<- paste("Est%",pc,sep="")
#GTL<-lntbc(ys,pc,cl)[1]; pGTL<-paste("UTL",pc,cl,sep="")
if( mth==1 ){GTL<-xtmp[3];pGTL<-paste("UTLa",pc,cl,sep="") }
if( mth==2 ) {GTL<-xtmp[5]; pGTL<-paste("UTLe",pc,cl,sep="") }
zOEL<- (log(OEL)-ys[1,1])/ys[1,2]; noel<-paste("z_OEL-",OEL,sep="")
out<-c(ys[1,1],ys[2,1],ys[1,2],ys[2,2],GM,GSD,EX,EXL,EXU,
      pc.obs,pc.est,GTL,zOEL)
onames<- c("mu","se.mu","sigma","se.sigma","GM","GSD","EX",nEXL,nEXU,
      pco,pce,pGTL,noel)
xmax<-max(d[,1])
m5 <- nptl(nt,pc/100,cl/100)
if( is.na(m5) )NPTL<-NA    else NPTL<- rev(sort(d[,1]))[m5]
k1<-  max(d[,1])*1.1
km<-kmms(d,cl)
nKML<-paste("KLCL",pc,sep="-"); nKMU<-paste("KUCL",pc,sep="-")
nNPTL<-paste("NpUTL",pc,cl,sep="")
LNstat<- c(out,km[1],km[3],km[4],NPTL,xmax,round(PCndet,2),round(nt),Rsq,ndet)
names(LNstat) <- c(onames,
 "KMmean",nKML,nKMU,nNPTL,"Maximum","NonDet%","n","Rsq","m")
#LNstat<-cbind(LNstat)
LNstat
}
############################ qqlognA  ####################################

qqlognA<-function(d=ex1,Ip="NONE",OEL=3000,unit="MSv x 100",cl=95,mth=1,ro=3){
#  qqlognA 4 OCT 2003 Maximum likelihood Estimates and
#  q-q plot for Left censored sample from  Lognormal Distribution 
# INPUT:d two column matris with
#         nonnegative data in d[,1] 
#         censoring(0=non-detect;1=detect) indicator in d[,2]
#       Ip  part of title used in plots ( default="NONE")
#       OEL is Ocupational Exposure Limit in data units*
#       unit is data units--default micrograms/m^3 (mug/m^3 
#       cl is confidence level for one-sided confidence intervals
#       ro controls rounding ( digidt to right of decimal point)
# REQUIRES:  lnstats() censored data Lognormal and Non parametric analysis
#            AND all functions listed in lnstats()  
du<-d[,1:2]
#Ip<- paste(Ip, deparse(substitute(d) )," (",unit,")" )
if(Ip=="NONE") Ip<-paste("Lognomal Q-Q Plot for",deparse(substitute(d) ) )
Ip<- paste(Ip, " (",unit,")" )

par(mfrow = c(1, 1), oma = c(1.1, 1.1, 1.1, 1.1), cex = 1.2)
kmg<- plend(du)
ym <- kmg[,2]  
ss<-lnstats(du,OEL,pc=95,cl,mth)
ss<-as.list(round(ss,ro))
Rsq<-round(ss$Rsq,3)
xq<- qnorm(kmg$apl)
pma<-max(1.1*kmg[,2])
         #xlim=c(-3,3),ylim=c(0,pma),
  plot( xq,ym,type = "n",xlim = c(min(xq),max(xq)),
        xlab="Normal Quantile",ylab=paste(" ",unit),log="y" )
points(xq,ym, pch = 1, cex = 0.6,col="red")
xx<- c(xq[1],xq[length(xq)]);  yhat<- exp(ss$mu + xx*ss$sigma)
lines(cbind(xx,yhat),type="b")

# Nef = solve sem= SD/sqrt(Nef) for Nef

tfile<-paste(   
   "\n GM= ",round(ss$GM,ro),"GSD= ",round(ss$GSD,ro),
   "\nArithmetic Mean ( ",cl,"% UCL ) ",
   "\n ", round( ss$EX,ro) ,"  (", round(ss$LCL,ro),
   ",",round( ss$UCL ,ro)," )",
   "\nKaplan-Meier Mean (",cl,"% UCL) ",
   "\n", ss$KMmean,"  (",ss$KLCL,",", ss$KUCL,")",
    "\nNon-Detects Percent= ",ss$NonD ,
	"\nSample Size n=", ss$n," m=",ss$m )   
tadd<- paste("\nZ for OEL",OEL,unit," =  ", round( (log(OEL) -ss$mu)/ss$sigma,2)) 
if(OEL > 0) tfile<-paste(tfile,tadd)
 	xp<- min(xq) +  0.2* abs( max(xq) - min(xq))
	yp<- min(ym) + 0.4*( max(ym) - min(ym))
 text(xp,yp,tfile,cex=0.7)
	tfile2<-paste(
        "\nObs.95th.Percentile= ", round(ss[[10]],ro), 
	"\nEst.95th.Percentile= ",round( ss[[11]],ro),
	"\n95-",cl," Geometric TL= ", round(ss$UTL,ro), 
	"\n 95-",cl,"Nonparam  TL=  ", ss$NpUTL,
	"\nMaximum Value = ",max(d[,1]),
	"\nR Square ", Rsq )
         xp<- 0.7*max(xq)
	yp<- min(ym) + 0.05*(max(ym) - min(ym))
text(xp,yp,tfile2,cex=0.7)
	Ip<-paste(Ip,"\n   Lognormal(",ss$mu,",",ss$sigma,")",
         "  Q-Q Plot ML Method=",mth,"Confidence Limits" )
#	"MLEs Method",mth,"Confidence Limits" )

	mtext( side = 3, line = 0, Ip , cex = 0.9,col="darkgreen")
	t1<-paste("qqlognA ",substring(date(),1,10),substring(date(),20,24))
	#mtext( side = 1, line = 4, t1 , cex = 0.7)
        
#INPUT(d =snl ,Ip="TITLE" ,OEL=0,unit="data units",cl=95,ro)
#
}

############################ lnexh1 ##################################
mlnd2 <- function(dd =ex1){
#      Exhibit 1 ORNL/TM-2004/146  R function mlnd2
#    find  Maximum Likelihood estimates of mu and sig 
#    left censored sample from lognormal distribution
# 
# USAGE: mlnd2(dd=ex1)
# ARGUMENT: matrix dd with d[i] in column 1 and cen in col 2
#  d[i] is positive lognormal data cen=0 for non-detect ; 1 for detect
#     y= log(d) is normal with mean mu and standard deviation sig
# VALUE: ML estimate of mu  and sig  estimates of standard errors
#          of mu and sig and - 2*Log-likelihood function
#         
# REQUIRES: function nlnd() and base R optim()   
#  initial estimate of mu and sig required by optim()
   yt <- ifelse(dd[,2]==0,dd[,1]/2,dd[,1])
   est <- c( mean(log(yt)), sd(log(yt) ) )
   n <- dim(dd)[[1]] ; m <- sum(dd[,2])
# ML estimate mu and sig
#
est <- optim(est,nlnd, method = c("Nelder-Mead"),dx=dd )$par
cont <- list(parscale=abs(est))
opt1 <- optim(est,nlnd ,NULL,  method ="L-BFGS-B",lower=c(-Inf,0.0),
          upper=c(Inf,Inf),cont, hessian=T,dx=dd )
mle <- opt1$par               #  ML estimates of mu and sig
vcm <- solve(opt1$hessian)    #  ML varaince-covariance matrix
semle <- sqrt(diag(vcm))      #  standard errors of mu and sig
cv <- vcm[1,2]                #  covariace(mu,sig) 
est <- c( round( c(mle,semle,cv) ,5 ),n,2*opt1$value)
names(est) <- c("mu","sig","se.mu","se.sig","cov","n","-2Log(L)")
est
}

nlnd <- function(p=est,dx){
# compute  - log liklihood for lognormal sample
 mu <- p[1]; sig <- p[2]; d <- dx[,1]
xx<-ifelse(dx[,2]==1,dlnorm(d,mu,sig,log=T) , plnorm(d,mu,sig,log.p=T))
  -sum(xx)
}

############################ lnexh2   ########################## 

lnexh2 <- function(ww=ex3,lod=30) {
#     Exhibit 2  ORNL/TM-2004/146
#    find ML estimates 1956-65 quarterly data
#   data in ww Col 1   Col 2    Col 3 
#              dose   cen(0,1)  t61
#     y = log(dose)  x = year - 1961
#       E(y) = alpha + beta*x
#
#    initial estimates using LS with zeros = lod/2
#
y0 <- log(ifelse( ww[,2]==0,lod/2,ww[,1] ))
go<-  summary( lm( y0  ~ ww[,3]  ) )
est <- c(go$coef[1,1], go$coef[2,1] , go$sigma ) 
names(est) <- c("alpha","beta","sigma")
#   Use R function optim() with lognomal log-likelihood LNlr2()
#   use Nelder-mead option to obtain ML estimates
opt <- optim(est,LNlr2, method = c("Nelder-Mead"),w=ww)
est <- opt$par
#   use "L-BFGS-B" option to obtain estimate of var-covar matrix
opt <- optim(est,LNlr2 ,NULL,  method ="L-BFGS-B",
    lower=c(-100.0,-100.0,0.1), upper=c(Inf,Inf,Inf),,hessian=T,w=ww)
# use "L-BFGS-B" option again with pscale added
pscale <- c(sqrt(diag(solve(opt$hessian))))
opt <- optim(est,LNlr2 ,NULL,  method ="L-BFGS-B",
    lower=c(-100.0,-100.0,0.1), upper=c(Inf,Inf,Inf),
    control=list(parscale=pscale), hessian=T,w=ww)
est <- opt$par
se <- sqrt(diag(solve(opt$hessian)))
vcvmle <- solve(opt$hessian) # ML covariance matrix
drc <- diag( 1/se) ; 
corm <- round( drc%*%vcvmle%*%drc,6) # correlation matrix
vcvmle <- round(vcvmle,9)

est <- c(opt$par,2*opt$value,opt$conver)
names(est) <- c("alpha","beta","sigma","-2Log(L)","Convrg")
vcv.cor <- cbind(vcvmle,corm)
se <- c(se,length(ww[,1]),NA)
mle <- rbind( est, se)
out <- list( mle ,vcvmle )
names(out) <- list("MLE","Variance-Covariance")
out
}

LNlr2<-function(par=est,w){
#  LNlr2 =  - log liklihood for left censored  sample lognormal
#
d<-w[,1]; cen<-w[,2] ; m<-par[1] + par[2]*w[,3] ;s<-par[3]
xx<-ifelse( cen==1,dlnorm(d,m,s,log=T) ,plnorm(d,m,s,log=T) )
  -sum(xx )
}

############################## readss.R ##########################

readss<-function(fn="Ex1") {
# to read fn.txt into R
# fn.txt may be space delimited file from Excel
#   with valid R column names
# USAGE: readss("fn")
# VALUE: data.frame 
#        writes output from lnstats to file "fnout.txt"
# 
# 
tmp<- read.table( paste(fn,"txt",sep=".") ,T )
stats<-lnstats(tmp,OEL=1250)
stats<-cbind(stats)
nout<-paste(fn,"out.txt",sep="")
sink(nout)
print(stats)
sink()
tmp
}

########################## end util.R ##########################