Nothing
mlefit<-function(x, dist="weibull", debias=NULL, optcontrol=NULL) {
## tz is required for MLEloglike and MLEsimplex calls now
default_tz=0
## sign is now required for MLEloglike call
default_sign=1
## check basic parameters of x
if(class(x)!="data.frame") {stop("mlefit takes a structured dataframe input, use mleframe")}
if(ncol(x)!=3) {stop("mlefit takes a structured dataframe input, use mleframe")}
xnames<-names(x)
if(xnames[1]!="left" || xnames[2]!="right"||xnames[3]!="qty") {
stop("mlefit takes a structured dataframe input, use mleframe") }
## test for any na's and stop, else testint below will be wrong
## initialize counts at zero, to be filled as found
Nf=0
Ns=0
Nd=0
Ni=0
## need this length information regardless of input object formation
failNDX<-which(x$right==x$left)
suspNDX<-which(x$right<0)
Nf_rows<-length(failNDX)
if(Nf_rows>0) {
Nf<-sum(x[failNDX,3])
}
Ns_rows<-length(suspNDX)
if(Ns_rows>0) {
Ns<-sum(x[suspNDX,3])
}
discoveryNDX<-which(x$left==0)
Nd_rows<-length(discoveryNDX)
if(Nd_rows>0) {
Nd<-sum(x[discoveryNDX,3])
}
testint<-x$right-x$left
intervalNDX<-which(testint>0)
interval<-x[intervalNDX,]
intervalsNDX<-which(interval$left>0)
Ni_rows<-length(intervalsNDX)
if(Ni_rows>0) {
Ni<-sum(x[intervalsNDX,3])
}
## further validate the input arguments for non-fsiq object
## somehow this code wrecked function on frames attributed with fsiq from mleframe
## if(length(attributes(x)$fsiq)!=1) {
## stop if Nf+Ns+Ndi != nrow(x)
## if( (Nf+Ns+Nd+Ni) != nrow(x)) {
## stop("invalid input dataframe")
## }
## rebuild input vector from components, just to be sure
## fsiq<-rbind(x[failNDX,], x[suspNDX,], x[discoveryNDX,], interval[intervalsNDX,])
## end input validation code
## }else{
## fsiq<-x
## }
## rebuild input vector from components, because this order is critical
fsiq<-rbind(x[failNDX,], x[suspNDX,], x[discoveryNDX,], interval[intervalsNDX,])
## Not sure what to place as restriction for C++ call
## if((Nf+Ni)<3) {stop("insufficient failure data")}
## now form the arguments for C++ call
## fsdi is the time vector to pass into C++
## data_est is used to estimate the magnitude of data
fsd<-NULL
data_est<-NULL
if((Nf+Ns)>0) {
fsd<-fsiq$left[1:(Nf_rows + Ns_rows)]
## assure that data_est is a clone
data_est<-fsiq$left[1:(Nf_rows + Ns_rows)]
}
if(Nd>0) {
fsd<-c(fsd,fsiq$right[(Nf_rows + Ns_rows + 1):(Nf_rows + Ns_rows + Nd_rows)])
data_est <- c(data_est, 0.5*(fsiq$right[(Nf_rows + Ns_rows + 1):(Nf_rows + Ns_rows + Nd_rows)]))
}
if(Ni>0) {
fsdi<-c(fsd, fsiq$left[(Nf_rows + Ns_rows + Nd_rows + 1):nrow(fsiq)],
fsiq$right[(Nf_rows + Ns_rows + Nd_rows + 1):nrow(fsiq)])
data_est<-c(data_est, (fsiq$left[(Nf_rows + Ns_rows + Nd_rows + 1):nrow(fsiq)] +
fsiq$right[(Nf_rows + Ns_rows + Nd_rows + 1):nrow(fsiq)])/2)
}else{
fsdi<-fsd
data_est<-fsd
}
q<-fsiq$qty
## third argument will be c(Nf,Ns,Nd,Ni)
N<-c(Nf_rows,Ns_rows,Nd_rows,Ni_rows)
## establish distribution number
if(tolower(dist)=="weibull" || tolower(dist)=="weibull2p" || tolower(dist)=="weibull3p") {
dist_num=1
m <- mean(log(data_est))
v <- var(log(data_est))
shape <- 1.2/sqrt(v)
scale <- exp(m + 0.572/shape)
vstart <- c(shape, scale)
}else{
if(tolower(dist)=="lognormal"|| tolower(dist)=="lognormal2p"|| tolower(dist)=="lognormal3p") {
dist_num=2
ml <- mean(log(data_est))
sdl<- sd(log(data_est))
vstart<-c(ml,sdl)
}else{
stop("distribution not resolved")
}
}
## Optional optimization control list to be handled here
## vstart will be as estimated
limit<-1e-5
maxit<-100
listout<-FALSE
if(length(optcontrol)>0) {
if(length(optcontrol$vstart>0)) {
vstart<-optcontrol$vstart
}
if(length(optcontrol$limit)>0) {
limit<-optcontrol$limit
}
if(length(optcontrol$maxit)>0) {
maxit<-optcontrol$maxit
}
if(length(optcontrol$listout)>0) {
listout<-optcontrol$listout
}
}
pos<-1
Q<-sum(q)
for(j in seq(1,4)) {
if(N[j]>0) {
Q<-c(Q, sum(q[pos:(pos+N[j]-1)]))
pos<-pos+N[j]
}else{
Q<-c(Q, 0)
}
}
names(Q)<-c("n","fo", "s", "d", "i")
MLEclassList<-list(fsdi=fsdi,q=q,N=N)
## Test for successful log-likelihood calculation with given vstart
## tz is required for MLEloglike call now
LLtest<-.Call("MLEloglike",MLEclassList,vstart,dist_num, default_sign, default_tz, package="abremDebias")
if(!is.finite(LLtest)) {
stop("Cannot start optimization with given parameters")
}
ControlList<-list(dist_num=dist_num,limit=limit,maxit=maxit)
## here is a good place to validate any debias argument (before more calculations begin)
if(length(debias)>0 && dist_num==1) {
if(tolower(debias)!="rba"&&tolower(debias)!="mean"&&tolower(debias)!="hirose-ross") {
stop("debias method not resolved")
}
}
## Handle the original 2 parameter case first
if(tolower(dist)=="weibull" || tolower(dist)=="lognormal" ||tolower(dist)=="weibull2p" || tolower(dist)=="lognormal2p" ) {
## listout control is passed as an integer to C++, this enables temporary change of status without losing input argument value
if(listout==TRUE) {
listout_int<-1
}else{
listout_int<-0
}
## tz inserted here with a default of zero
result_of_simplex_call<-.Call("MLEsimplex",MLEclassList, ControlList, vstart, default_tz, listout_int, package="abremDebias")
## extract fit vector from result of call to enable finishing treatment of the outvec
if(listout==FALSE) {
resultvec<-result_of_simplex_call
}else{
resultvec<-result_of_simplex_call[[1]]
}
outvec<-resultvec[1:3]
if(resultvec[4]>0) {
warn<-"likelihood optimization did not converge"
attr(outvec,"warning")<-warn
}
if(dist_num == 1) {
names(outvec)<-c("Eta","Beta","LL")
if(length(debias)>0) {
if(debias!="rba"&&debias!="mean"&&debias!="hirose-ross") {
stop("debias method not resolved")
}
if(debias=="rba") {
outvec[2]<-outvec[2]*rba(Q[1]-Q[3], dist="weibull",basis="median")
}
if(debias=="mean") {
outvec[2]<-outvec[2]*rba(Q[1]-Q[3], dist="weibull",basis="mean")
}
if(debias=="hirose-ross") {
outvec[2]<-outvec[2]*hrbu(Q[1]-Q[3], Q[3])
}
outvec[3]<-.Call("MLEloglike",MLEclassList,c(outvec[2],outvec[1]),dist_num, default_sign, default_tz, package="abremDebias")
attr(outvec,"bias_adj")<-debias
}
}
if(dist_num == 2) {
names(outvec)<-c("Mulog","Sigmalog","LL")
if(length(debias)>0) {
outvec[2]<-outvec[2]*rba(Q[1]-Q[3], dist="lognormal")
if(debias!="rba") {
warning("rba has been applied to adjust lognormal")
debias="rba"
}
outvec[3]<-.Call("MLEloglike",MLEclassList,c(outvec[1],outvec[2]),dist_num, default_sign, default_tz, package="abremDebias")
attr(outvec,"bias_adj")<-debias
}
}
if(listout==TRUE) {
optDF<-as.data.frame(result_of_simplex_call[[2]])
if(dist_num == 1) {
names(optDF)<-c("beta_est", "eta_est", "negLL", "error")
}
if(dist_num == 2) {
names(optDF)<-c("mulog_est", "sigmalog_est", "negLL", "error")
}
}
## end of 2p code
}
## this section of code is specifically addressing 3p models
if(tolower(dist)=="weibull3p" || tolower(dist)=="lognormal3p" ) {
## For now, listout is passed as integer
## listout argument has different meaning for 3p models
listout_int<-0
## for now enter a default tz=0
result_of_simplex_call<-.Call("MLEsimplex",MLEclassList, ControlList, vstart, default_tz, listout_int, package="abremDebias")
if(result_of_simplex_call[4]>0) {
stop("2p model does not converge")
}
## restore the meaning of listout
if(listout==TRUE) {
listout_int<-1
}else{
listout_int<-0
}
##########################################################
## this would be the position to go to C++
##########################################################
secant_warning=FALSE
## Tao Pang's original variable labels from FORTRAN are used where possible
DL<-ControlList$limit
## Introduce constraints for tz
## need to create the vector of potential minimums
fdr<-NULL
if(Nf>0) {fdr<-fsdi[1:Nf]}
if(Nd>0) {fdr<-c(fdr,fsdi[(Nf+Ns+1):(Nf+Ns+Nd)])}
if(Ni>0) {fdr<-c(fdr, fsdi[(Nf+Ns+Nd+Ni+1):(Nf+Ns+Nd+2*Ni)])}
C1<-min(fdr)
maxit<-100
## initial step is based on min(x)*.1
DX<-C1*0.1
X0<-0.0
istep<-0
X1<-X0+DX
if(X1>C1) {X1<-X0+0.9*(C1-X0)}
tz=0
FX0vec<-.Call("MLEdMaxLLdx", MLEclassList, ControlList, vstart, tz, package="abremDebias")
FX0<-FX0vec[1]
## new start estimate from last fit (without any modification)
vstart<-FX0vec[-1]
## X1 is next proposed tz
tz=X1
FX1vec<-.Call("MLEdMaxLLdx", MLEclassList, ControlList, vstart, tz, package="abremDebias")
FX1<-FX1vec[1]
## new start estimate from last fit (without any modification)
vstart<-FX1vec[-1]
## FX1 will contain slope sign information to be used only one time to find X2
D<- abs(FX1-FX0)
X2<-X1+abs(X1-X0)*FX1/D
if(X2>C1) {X2<-X1+0.9*(C1-X1)}
X0<-X1
X1<-X2
DX<-X1-X0
istep<-istep+1
## Detail output to be available with listout==TRUE
DF<-data.frame(steps=istep,root=X0,error=DX,deriv=FX1)
while(abs(DX)>DL&&istep<maxit) {
FX0<-FX1
## save last successful call with useful vstart
FX0vec<-FX1vec
## X1 is next proposed tz
tz=X1
FX1vec<-.Call("MLEdMaxLLdx", MLEclassList, ControlList, vstart, tz, package="abremDebias")
FX1<-FX1vec[1]
if(is.nan(FX1)) {
FX1<-FX0
secant_warning=TRUE
break
}
## new start estimate from last fit (without any modification)
vstart<-FX1vec[-1]
## FX1 will contain slope information only one time
D<- abs(FX1-FX0)
X2<-X1+abs(X1-X0)*FX1/D
if(X2>=C1) {X2<-X1+0.9*(C1-X1)}
X0<-X1
X1<-X2
DX<-X1-X0
istep<-istep+1
DFline<-data.frame(steps=istep,root=X0,error=DX,deriv=FX1)
DF<-rbind(DF,DFline)
## return to while loop
}
## provide a last good vstart in case FX1vec indicated nan
vstart<-FX0vec[-1]
## Can X0 be first trial, but ultimately subject to convergence problems??
listout_int<-0
result_of_simplex_call<-.Call("MLEsimplex",MLEclassList, ControlList, vstart, X0, listout_int, package="abremDebias")
## extract fit vector from result of call to enable finishing treatment of the outvec
outvec<-c(result_of_simplex_call[1:2], X0, result_of_simplex_call[3])
if(dist_num==1) {
names(outvec)<-c("Eta","Beta", "t0", "LL")
if(length(debias)>0) {
if(debias=="rba") {
outvec[2]<-outvec[2]*rba(Q[1]-Q[3], dist="weibull",basis="median")
}
if(debias=="mean") {
outvec[2]<-outvec[2]*rba(Q[1]-Q[3], dist="weibull",basis="mean")
}
if(debias=="hirose-ross") {
outvec[2]<-outvec[2]*hrbu(Q[1]-Q[3], Q[3])
}
outvec[3]<-.Call("MLEloglike",MLEclassList,c(outvec[2],outvec[1]),dist_num, default_sign, X0, package="abremDebias")
attr(outvec,"bias_adj")<-debias
}
}
if(dist_num == 2) {
names(outvec)<-c("Mulog","Sigmalog", "t0", "LL")
if(length(debias)>0) {
outvec[2]<-outvec[2]*rba(Q[1]-Q[3], dist="lognormal")
if(debias!="rba") {
warning("rba has been applied to adjust lognormal")
debias="rba"
}
outvec[3]<-.Call("MLEloglike",MLEclassList,c(outvec[1],outvec[2]),dist_num, default_sign, X0, package="abremDebias")
attr(outvec,"bias_adj")<-debias
}
}
optDF<-DF[-1]
## end of 3p code
}
## the following applies to both 2p and 3p results
attr(outvec,"data_types")<-Q[-2]
if(listout==FALSE) {
out_object<-outvec
}else{
out_object<-list(fit=outvec, opt=optDF)
}
out_object
## end function
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.