In Auburngrads/SMRD: Statistical Methods For Reliability Data

SMRD:::vinny(cache = T)
library(SMRD)

In this echapter...

Warning: some of the ML estimation attempts fail to converge

LocomotiveControl.ld <- frame.to.ld(locomotivecontrol,
                                    response.column = 1, 
                                    censor.column = 2,
                                    case.weight.column = 3,
                                    time.units = "kMiles")

LocomotiveControl.gamma.gmle.out <- Gamma.mle(LocomotiveControl.ld)

gmleprobplot(LocomotiveControl.ld,
             distribution = "gamma", 
             xlim = c(10,199),
             ylim = c(.0011,.991))

These examples can take longer to run

## gamma dist examples
lzbearing.ld <- frame.to.ld(lzbearing, response.column = 1)

bear.gamma.gmle.out <- Gamma.mle(lzbearing.ld)

summary(bear.gamma.gmle.out)

gmleprobplot(lzbearing.ld,
             distribution = "gamma",
             xlim = c(10,199),
             ylim = c(.0011,.991),
             compare = c("Lognormal","Weibull"))

legend(x = c(3.6,5.2),
       y = c(-1.6,-2.7), 
       c("Gamma","Lognormal","Weibull","95% pointwise confidence intervals"),
       lty = c(1,3,4,2),
       bty = "n",
       cex = 1.2)

bkfatigue10.ld <- frame.to.ld(bkfatigue10, 
                              response.column = 1,
                              time.units = "Kilocycles")

summary(bkfatigue10.ld)

bkfatigue10.gamma.gmle.out <- Gamma.mle(bkfatigue10.ld)

gmleprobplot(bkfatigue10.ld,
             distribution = "gamma",
             compare = c("Lognormal"))

gmleprobplot(bkfatigue10.ld,
             distribution = "Lognormal",
             compare = c("gamma"))

at7987.ld <- frame.to.ld(at7987, 
                         response.column = 1,
                         censor.column = 2, 
                         case.weight.column = 3,
                         time.units = "Kilocycles")

summary(at7987.ld)

at7987.gamma.gmle.out <- Gamma.mle(at7987.ld)

gmleprobplot(at7987.ld,
             distribution = "Gamma",
             compare = c("Lognormal"))

egeng examples

at7987.egeng.gmle.out <- egeng.mle(at7987.ld)

bear.egeng.gmle.out <- egeng.mle(lzbearing.ld)

lzbearing5.ld <- small.interval(lzbearing.ld, delta = 1e-03, tlog = T)

bear5.egeng.gmle.out <- egeng.mle(lzbearing5.ld)

at7987.egeng.gmle.out <- gmleprobplot(at7987.ld,
                                      distribution = "egeng")

gmleprobplot(lzbearing.ld,
             distribution = "egeng",
             xlim = c(10,199),
             ylim = c(.0011,.991),
             compare = c("Weibull","Lognormal","Exponential"))

text(2.4969, -1.31908, "EXP")
text(2.50633, -3.4758, "WEIB")
text(2.34606, -5.75833, "EGENG")
text(2.70432, -6.94453, "LOGNOR")

tmp <- one.dim.profile(bear.egeng.gmle.out,
                       profile.on.list = 3,
                       range.list = list(c(-1,3.0)),
                       size = 50,
                       save.structures = F)

profile.plot(tmp, variable.name = "lambda")
abline(v = 0, lty = 3)
abline(v = 1, lty = 3)
text(0.829561, 0.224463, "WEIB")
text(-0.236364, 0.224463, "LOGNOR")

Fan.ld <- frame.to.ld(fan,
                      response.column = 1, 
                      censor.column = 2, 
                      case.weight.column = 3,
                      time.units = "Hours")

Fan.egeng.gmle.out <- egeng.mle(Fan.ld)

gmleprobplot(Fan.ld,
             distribution = "egeng",
             xlim = c(200,9999),
             ylim = c(.0011,.39),
             compare = c("Weibull","Lognormal","Exponential"))

text(5.37368, -2.28317, "EXP")
text(5.39218, -2.6004, "WEIB")
text(6.00877, -3.03298, "GENG")
text(5.52783, -2.91186, "LOGNOR")


one.dim.profile(Fan.egeng.gmle.out,
                profile.on.list = 3,
                range.list = list(c(-8,2)),
                size = 50)

# need to rename this function to s3 method plot.one.dim.out
# This is interfering with s3 method stats::profile
SMRD:::profile.plot(Fan.egeng.gmle.outstruct3)

bisa examples

bear.bisa.gmle.out <- bisa.mle(lzbearing.ld)

gmleprobplot(lzbearing.ld,
             distribution = "bisa",
             compare = c("Lognormal"))

Fan.bisa.gmle.out <- bisa.mle(Fan.ld)

basic.gmleprobplot(Fan.ld,
                   distribution = "bisa",
                   xxx.mle.out = Fan.bisa.gmle.out,
                   compare = c("Lognormal"))

at7987.bisa.gmle.out <- bisa.mle(at7987.ld)

basic.gmleprobplot(at7987.ld,
                   distribution = "bisa",
                   xxx.mle.out = at7987.bisa.gmle.out,
                   compare = c("Lognormal"))

basic.gmleprobplot(at7987.ld,
                   distribution = "bisa",
                   plot.dist = "lognormal",
                   xxx.mle.out = at7987.bisa.gmle.out,
                   compare = c("Lognormal"))

igau examples

bear.igau.gmle.out <- igau.mle(lzbearing.ld)

bkfatigue10.igau.gmle.out <- igau.mle(bkfatigue10.ld)

basic.gmleprobplot(bkfatigue10.ld, 
                   distribution = "igau",
                   xxx.mle.out = bkfatigue10.igau.gmle.out,
                   compare = c("Lognormal"))

gmleprobplot(bkfatigue10.ld, 
             distribution = "igau",
             compare = c("Lognormal"))

gmleprobplot(bkfatigue10.ld,
             distribution = "bisa",
             compare = c("Lognormal"))
gmleprobplot(bkfatigue10.ld, 
             distribution = "gamma",
             compare = c("Lognormal"))

gmleprobplot(bkfatigue10.ld,
             distribution = "lognormal",
             compare = c("Lognormal"))

gmleprobplot(bkfatigue10.ld, 
             distribution = "Weibull",
             compare = c("Lognormal"))

gets (threshold) examples

bear.gets.sev.gmle.out <- gets.mle(lzbearing.ld, distribution = "sev")

Fan.gets.sev.gmle.out <- gets.mle(small.interval(Fan.ld, delta = 0.01),
                                  distribution = "sev")

Fan.gets.nor.gmle.out <- gets.mle(small.interval(Fan.ld, delta = 0.01),
                                  distribution = "normal")

basic.gmleprobplot(small.interval(Fan.ld, delta = 0.01),
                   distribution = "sevgets",
                   xxx.mle.out = Fan.gets.sev.gmle.out, 
                   my.title = "", 
                   cexlab = 1.5,
                   compare = list(Fan.gets.nor.gmle.out))

basic.gmleprobplot(small.interval(Fan.ld, delta = .01), 
                   distribution = "sevgets",
                   xxx.mle.out = Fan.gets.sev.gmle.out, 
                   my.title="", 
                   cexlab = 1.5,
                   compare = c("lognormal","gamma","bisa"))

gmleprobplot(small.interval(Fan.ld, delta = 0.01), 
             distribution = "sevgets",
             my.title = "", 
             cexlab = 1.5)

at7987.sevgets.gmle.out <- gets.mle(at7987.ld, distribution = "sev")

at7987.norgets.gmle.out <- gets.mle(at7987.ld, distribution = "normal")

basic.gmleprobplot(at7987.ld,
                   distribution = "sevgets",
                   xxx.mle.out = at7987.sevgets.gmle.out, 
                   my.title = "",
                   cexlab = 1.5,
                   compare = list(at7987.norgets.gmle.out))

# problem with data
AlloyC.ld <- frame.to.ld(alloyc,
                         response.column = c(1,2),
                         censor.column = 3,
                         case.weight.column = 4,
                         data.title = "Alloy C",
                         time.units = "ksi")

AlloyC.sevgets.gmle.out <- gets.mle(AlloyC.ld, distribution = "sev")

AlloyC.norgets.gmle.out <- gets.mle(AlloyC.ld, distribution = "normal")

basic.gmleprobplot(AlloyC.ld, distribution = "sevgets",
                   xxx.mle.out = AlloyC.sevgets.gmle.out,
                   my.title = "",
                   cexlab = 1.5,
                   plot.dist = "sev")

AlloyC.gets.nor.gmle.out <- gets.mle(AlloyC.ld, distribution = "normal")

basic.gmleprobplot(AlloyC.ld,
                   distribution = "normalgets",
                   xxx.mle.out = AlloyC.norgets.gmle.out, 
                   my.title = "", 
                   cexlab = 1.5,
                   plot.dist = "sev")

Begin to fill region for lr confidence intervals. The code in the following three chunks is experimental

Fan.egeng.gmle.out <- FillRegion(Fan.egeng.gmle.out,
                                 nbound = 10,
                                 iter = 500)

summary(Fan.egeng.gmle.out.jcr)

names(Fan.egeng.gmle.out.jcr)

basic.gmleprobplot(Fan.ld,distribution = "egeng",
                   xlim = c(200,99999),
                   ylim = c(.0011,.69),
                   xxx.mle.out = Fan.egeng.gmle.out,
                   my.title = "",
                   cexlab = 1.5,
                   conlev = .95,
                   ciMethod = "lr.approx",
                   length.time.vec = 2)

bear.egeng.gmle.out <- egeng.mle(lzbearing.ld)


Fan.weibull.gmle.out <- ls2.mle(Fan.ld, distribution = "weibull")


Fan.weibull.gmle.out <-  FillRegion(Fan.weibull.gmle.out,
                                    nbound = 4,
                                    iter = 50,
                                    cull = 2 )

summary(Fan.weibull.gmle.out.jcr)
summary(Fan.weibull.gmle.out)

mleprobplot(Fan.ld,
            distribution = "Weibull",
            xlim = c(200,9999),
            ylim = c(.0031,.49))

basic.gmleprobplot(Fan.ld,
                   distribution = "Weibull",
                   xlim = c(200,9999),
                   ylim = c(.0031,.49),
                   xxx.mle.out = Fan.weibull.gmle.out,
                   my.title = "",
                   cexlab = 1.5,
                   conlev = .95,
                   length.time.vec = 30)

basic.gmleprobplot(Fan.ld,
                   distribution = "Weibull",
                   xlim = c(200,9999),
                   ylim = c(.0031,.49),
                   xxx.mle.out = Fan.weibull.gmle.out,
                   my.title = "",
                   cexlab = 1.5,
                   conlev = 0.95,
                   ciMethod = "lr.approx",
                   length.time.vec = 30)

Issues with this code chunk.

fcn = function(theta,time,distribution){

      SMRD:::wqmf.phibf((log(time)-theta[1]) / exp(theta[2]),
                      distribution = "Weibull")

}

Fr.conf(fcn, 
        fcn.arg2 = log.seq(200,2000,length=10),
        gmle.out = Fan.weibull.gmle.out,
        ptwise = T)

Fr.conf(fcn,
        fcn.arg2 = log.seq(200,2000,length = 10),
        gmle.out = Fan.weibull.gmle.out,  
        ptwise = T, 
        extrapolate = T)

Fan.weibull.gmle.out <- ls.mle(Fan.ld,
                               distribution = "Weibull")

Fan.lev.gmle.out <- ls.mle(Fan.ld,
                           distribution = "lev")
names(Fan.weibull.gmle.out)

 Fan.weibull.gmle.out <- FillRegion(Fan.weibull.gmle.out,
                                    nbound = 4,
                                    iter = 500, 
                                    cull = 2)

summary(Fan.weibull.gmle.out.jcr)

mleprobplot(Fan.ld, 
            distribution = "Lognormal")

basic.gmleprobplot(Fan.ld,
                   distribution = "Weibull",
                   plot.dist = "Weibull",
                   xxx.mle.out = Fan.weibull.gmle.out,
                   my.title = "",
                   cexlab = 1.5,
                   conf.level = 0.95,
                   ciMethod = "lr.approx",
                   length.time.vec = 20)

Fan.lognormal.gmle.out <- ls.mle(Fan.ld,
                                 distribution = "Lognormal")

basic.gmleprobplot(Fan.ld,
                   distribution = "Lognormal", 
                   plot.dist = "Lognormal",
                   xxx.mle.out = Fan.lognormal.gmle.out,
                   my.title ="",
                   cexlab = 1.5,
                   conlev = 0.95,
                   ciMethod = "lr.approx",
                   length.time.vec = 20)

BearingCage.ld <- frame.to.ld(bearingcage,
                              response.column = 1, 
                              censor.column = 2, 
                              case.weight.column = 3,
                              time.units = "Hours")

summary(BearingCage.ld)

BearingCage.weibull.gmle.out <- ls.mle(BearingCage.ld,
                                       distribution = "Weibull")

names(BearingCage.weibull.gmle.out)

BearingCage.weibull.gmle.out.jcr <- FillRegion(BearingCage.weibull.gmle.out,
                                               nbound = 4,
                                               iter = 500,
                                               cull = 2)

summary(BearingCage.weibull.gmle.out.jcr)

mleprobplot(BearingCage.ld,
            distribution = "Weibull",
            xlim = c(200,10000),
            ylim = c(.00031,.19),
            time.vec = log.seq(200,10000,20))

basic.gmleprobplot(BearingCage.ld,
                   distribution = "Weibull",
                   xxx.mle.out = BearingCage.weibull.gmle.out,
                   my.title = "",
                   cexlab = 1.5,
                   conf.level = 0.95,
                   length.time.vec = 20,
                   xlim = c(200,10000),
                   ylim = c(.00031,.19), 
                   time.vec = log.seq(200,10000,20),
                   ciMethod = "lr.approx")

basic.gmleprobplot(BearingCage.ld,
                   distribution = "Weibull",
                   xxx.mle.out = BearingCage.weibull.gmle.out,
                   my.title = "",
                   cexlab = 1.5,
                   conf.level = 0.95,
                   length.time.vec = 20, 
                   xlim = c(200,10000),
                   ylim = c(.00031,.19), 
                   time.vec = log.seq(200,10000,20))

SMRD:::profile.plot(Fr.profile(fcn <- function(x){x[1]},
                        BearingCage.weibull.gmle.out))

tmppro <- Fr.profile(fcn <- function(x){x[1]}, 
                     BearingCage.weibull.gmle.out)

#write a wrapper to specify the function for profile (e.g., orig param, param, #quantiles, failure prob)


vchan.functions <- c("Fr.profile",
                     "FillRegion",
                     "Fr.conf",
                     "summary.FillRegion.out",
                     "basic.gmleprobplot",
                     "gmleprobplot")

SMRD:::profile.plot(tmppro)

do.profile.plot(BearingCage.weibull.gmle.out,
                parameter = "log(scale)")

do.profile.plot(BearingCage.weibull.gmle.out,
                original.parameter = "scale")

do.profile.plot(BearingCage.weibull.gmle.out,
                quantile = 0.1)

do.profile.plot(BearingCage.weibull.gmle.out,
                failure.probability = 5000)

BearingCage.weibull.gmle.out$model$f.origparam

Right truncation example

doatrun.ld <- frame.to.ld(doatrun,
                          response.column = c(1,2),
                          censor.column = 3, 
                          case.weight.column = 4,
                          truncation.response.column = 5, 
                          truncation.type.column = 6, 
                          data.title = "DOA Truncated Data", 
                          time.units = "Hours")

summary(doatrun.ld)

cdfest(doatrun.ld)

mlest(doatrun.ld,"Weibull")

mleprobplot(doatrun.ld,
            distribution = "Weibull")

mleprobplot(doatrun.ld,
            distribution = "Weibull",
            trunc.correct = F)

npprobplot(doatrun.ld,
           distribution = "Weibull")

npprobplot(doatrun.ld,
           distribution = "Weibull",
           trunc.correct = F)

left truncation example

lfptrun100.ld <- frame.to.ld(lfptrun100, 
                             response.column = c(1), 
                             case.weight.column = 2,
                             truncation.response.column = 3, 
                             truncation.type.column = 4, 
                             data.title = "LFP Truncated Data", 
                             time.units = "Hours")

summary(lfptrun100.ld)

cdfest(lfptrun100.ld)

mlest(lfptrun100.ld,"Weibull")

mleprobplot(lfptrun100.ld,
            distribution = "Weibull",
            trunc.correct = F)

mleprobplot(lfptrun100.ld,
            distribution = "Weibull")

plot(lfptrun100.ld,
     distribution = "Weibull",
     trunc.correct = F)

plot(lfptrun100.ld,
     distribution = "Weibull")

Examples of fitting the Random Fatigue Limit Model

Laminate Panel Fatigue Data

LaminatePanel.ld <- frame.to.ld(laminatepanel,
                                response.column = "kilocycles",
                                time.units = "kilocycles",
                                data.title = "Laminate Panel Fatigue Data",
                                censor.column = "event", 
                                x.column = "mpa")

censored.data.plot(LaminatePanel.ld, 
                   x.axis = "log", 
                   y.axis = "log",
                   response.on.yaxis = F)

censored.data.plot(LaminatePanel.ld)

LaminatePanel.fit11.out <- rflm.mle(LaminatePanel.ld,
                                    cond.dist = "sev", 
                                    fl.dist = "sev")

plot(LaminatePanel.fit11.out,
     response.on.yaxis = F)

LaminatePanel.fit11f.out <- rflm.mle(LaminatePanel.ld,
                                     cond.dist = "sev", 
                                     fl.dist = "sev",
                                     debug = 1,
                                     fixed.param.list =   
                                     list(fixed.parameters = 3,
                                          fixed.parameter.values = 0.20))

plot(LaminatePanel.fit11f.out,
     response.on.yaxis = F)


LaminatePanel.fit12.out <- rflm.mle(LaminatePanel.ld,
                                    cond.dist = "sev", 
                                    fl.dist = "normal")

print(LaminatePanel.fit12.out)

LaminatePanel.fit12.out <- rflm.mle(LaminatePanel.ld, 
                                    cond.dist = "normal", 
                                    fl.dist = "sev")

print(LaminatePanel.fit12.out)

plot(LaminatePanel.fit12.out, 
     response.on.yaxis = F)

LaminatePanel.fit22.out <- rflm.mle(LaminatePanel.ld,
                                    cond.dist = "normal", 
                                    fl.dist = "normal")

print(LaminatePanel.fit22.out)

The code in the following chunk is for testing only

one.dim.profile(LaminatePanel.fit22.out,
                profile.on.list = 3,
                size = 50, 
                range.list = list(c(-35,0)),
                save.parameter.vectors = T)

one.dim.profile(LaminatePanel.fit22.out,
                profile.on.list = 3,
                size = 300, 
                range.list = list(c(-35,0)), 
                save.parameter.vectors = T)

SMRD:::profile.plot(LaminatePanel.fit22.outstruct3)

one.dim.profile(LaminatePanel.fit22.out,
                profile.on.list = 5,
                size = 30, 
                range.list = list(c(-9.21, -2.3)),
                save.parameter.vectors = T)

inconel superalloy data

inconelsub <- inconel[inconel$Strain < 0.007,]

inconel.ld <- frame.to.ld(inconel,
                          response.column = "cycles",
                          data.title = "Inconel Superalloy Fatigue Data with Stress < 7000",
                          censor.column = "event", 
                          x.column = "strain")

censored.data.plot(inconel.ld, 
                   y.axis = "log",
                   response.on.yaxis = F)

inconel.fit11 <- rflm.mle(inconel.ld, 
                          cond.dist = "sev",
                          fl.dist = "sev")
print(inconel.fit11)

plot(inconel.fit11)

inconel.fit12 <- rflm.mle(inconel.ld, 
                          cond.dist = "sev",
                          fl.dist = "normal")
print(inconel.fit12)

inconel.fit21 <- rflm.mle(inconel.ld, 
                          cond.dist = "normal", 
                          fl.dist = "sev")#
print(inconel.fit21)

inconel.fit22 <- rflm.mle(inconel.ld, 
                          cond.dist = "normal", 
                          fl.dist = "normal")
print(inconel.fit22)

plot(inconel.fit22)

titanium Ti64 fatigue data

titanium.ld <- frame.to.ld(titanium,
                           response.column = "kilocycles",
                           censor.column = "event",
                           data.title = "Ti64 Fatigue Failure Data",
                           x.column = "strain")

censored.data.plot(titanium.ld,
                   xlim = c(5.1,24),
                   y.axis = "log",
                   x.axis = "log",
                   response.on.yaxis = F)

titanium.fit11 <-rflm.mle(titanium.ld, 
                          cond.dist = "sev",
                          fl.dist = "sev")
print(titanium.fit11)
plot(titanium.fit11)

titanium.fit12 <-rflm.mle(titanium.ld, 
                          cond.dist = "sev",
                          fl.dist = "normal")

print(titanium.fit12)

titanium.fit21 <- rflm.mle(titanium.ld, 
                           cond.dist = "normal", 
                           fl.dist = "sev")
print(titanium.fit21)

titanium.fit22 <- rflm.mle(titanium.ld, 
                           cond.dist = "normal", 
                           fl.dist = "normal")
print(titanium.fit22)
plot(titanium.fit22)

plot(titanium.fit22, 
     response.on.yaxis = F,
     xlim = c(5.1,24))

Concrete Data

concrete[1:5,]

Concrete.ld <- frame.to.ld(concrete,
                           response.column = "kilocycles",
                           time.units = "Kilocycles",
                           data.title = "Concrete Fatigue Data",
                           xlabel = "Smax-Sfail",
                           x.column = "ratio")

censored.data.plot(Concrete.ld, 
                   x.axis = "log", 
                   y.axis = "log",
                   response.on.yaxis = F)

Concrete.fit11.out <- rflm.mle(Concrete.ld,
                               cond.dist = "sev",
                               fl.dist = "sev")
print(Concrete.fit11.out)

plot(Concrete.fit11.out,
     response.on.yaxis = F)

Concrete.fit12.out <- rflm.mle(Concrete.ld,
                               cond.dist = "sev",
                               fl.dist = "normal")
print(Concrete.fit12.out)

Concrete.fit21.out <- rflm.mle(Concrete.ld,
                               cond.dist = "normal",
                               fl.dist = "sev")
print(Concrete.fit21.out)


Concrete.fit22.out <- rflm.mle(Concrete.ld,
                               cond.dist = "normal",
                               fl.dist = "normal")

print(Concrete.fit22.out)

plot(Concrete.fit22.out,
     response.on.yaxis = F)

high-cycle fatigue

hcf.ld <- frame.to.ld(hcfdata,
                      response.column = "cycles",
                      data.title = "High Cycle Fatigue Data",
                      censor.column = "event", 
                      x.column = c("swt"))

censored.data.plot(hcf.ld, 
                   x.axis = "log",
                   y.axis = "log")

censored.data.plot(hcf.ld, 
                   x.axis = "log", 
                   y.axis = "log", 
                   response.on.yaxis = F)

# Turning the trace on and up the number of digits for these

hcf.fit22f1 <- rflm.mle(hcf.ld, 
                        cond.dist = "normal", 
                        fl.dist = "normal",
                        digits = 10,
                        dump = 7,
                        fixed.param.list=list(fixed.parameters = 3,
                                              fixed.parameter.values = 0.02))
#likelihood on Alpha chip= -123.3193856
#likelihood on Wintel chip= -123.31939

hcf.fit22f2 = rflm.mle(hcf.ld, 
                       cond.dist = "normal",
                       fl.dist = "normal",
                       digits = 10,
                       dump = 7,
                       fixed.param.list = list(fixed.parameters = 3,
                                               fixed.parameter.values = 0.0002))
#likelihood on Alpha chip = -123.2937926
#likelihood on Wintel chip= -123.29379

hcf.fit22f3 = rflm.mle(hcf.ld, 
                       cond.dist = "normal", 
                       fl.dist = "sev",
                       digits = 10,
                       dump = 7,
                       fixed.param.list = list(fixed.parameters = 3,
                                               fixed.parameter.values = 0.00002))
# likelihood on Alpha chip= -123.2937897 (false convergence messaage,
# but the gradiant looks fine)
# likelihood on Wintel chip= -123.29620
#
# Notice that as the fixed value of sigma approaches zero, the likelihood
# in the constrained optimization stops increasing.
#
# This indicates a ridge in the likelihood surface.
# Using the last estimates as starting values for an unconstrained optimization
# does not yield any improvement, but we get a successful convergence
# (perhaps we should not) message.
#
# The following converges on my Alpha (128 bit double precision arithmetic)
# but not on Windows (64 bit double precision arithmetic)
#
# The Alpha gave a standard error of 0.02538 for an estimate of 2e-05!
# But no standard error from windows.
hcf22.fit <- rflm.mle(hcf.ld, 
                      cond.dist = "normal", 
                      fl.dist = "normal",
                      digits = 10,
                      dump = 7, 
                      theta.start = c(9.356277116, 
                                      -2.935035745, 
                                      2e-04, 
                                      3.863027231, 
                                      0.1183292607))

#likelihood on Alpha chip = -123.2937897
#likelihood on Wintel chip= -123.2937926

plot(hcf.fit22f3, 
     response.on.yaxis = F)

plot(hcf22.fit, response.on.yaxis = F)

# This take a while to run, but the result is interesting
#
#Warning: with the default size=40, the following command (uses a loop)
#crashed my windows session as I ran out of RAM (128MB).
#No problem on my 1GB Alpha. My guess is that it will run ok with 256MB
#It ran ok (about 15 minutes on my 233MHz) with  size = 20
censored.data.plot(hcf.ld, 
                   x.axis = "log", 
                   y.axis = "log", 
                   response.on.yaxis = F)

hcf.fit <- rflm.mle(hcf.ld, 
                    cond.dist = "sev",
                    fl.dist = "normal")

plot(hcf.fit)

# The following indicates that we will have trouble  
# identifying sigma from these data
one.dim.profile(hcf22.fit,
                profile.on.list = 3,
                size = 20, 
                range.list = list(c(-10,0)))

SMRD:::profile.plot(hcf22.fitstruct3)

#check
xcdf(1, 2, 29, -5, 270 ,  .36, 5.3, .020,6.0,1,0.0,0 )
#jave[1] 0.009247504
#>
xpdf(1, 2, 29, -5, 380 ,  .36, 5.3, .020,3.5322 )
#jave[1] 0.2656861

Auburngrads/SMRD documentation built on Sept. 14, 2020, 2:21 a.m.