Nothing
"mps2par" <-
function(x, type, init.para=NULL, ties=c("bernstein", "rounding", "density"),
delta=0, log10offset=3, get.untied=FALSE, check.support=TRUE,
moran=TRUE, silent=TRUE, null.on.not.converge=TRUE,
ptransf= function(t) return(t),
pretransf=function(t) return(t),
mle2par=TRUE, ...) {
x <- sort(x); n <- length(x); untied <- FALSE; ITS <- NA
ties.method <- match.arg(ties)
if(is.null(type)) {
warning("must specify distribution type")
return(NULL)
}
if(delta < 0) {
warning("a delta argument < 0 does not make sense")
return(NULL)
}
if(log10offset <= 0) {
warning("a log10offset argument <= 0 does not make sense")
return(NULL)
}
cfunc <- "1" # The following is experimental and following
# a paper in the literature: Ghosh, K., and Jammalamadaka, S.R., 2001,
# A general estimation method using spacings: Journal of Statistical Planning
# and Inference, v. 93, pp. 71--82.
the.cfunc <- switch(cfunc, "1" =function(x) -log(x) ,
"2" =function(x) log(-log(x)),
"3" =function(x) x*log(x) ,
"4" =function(x) x^2,
"5" =function(x) -sqrt(x) ,
"6" =function(x) 1/x ,
"7" =function(x) abs(x-1))
# The Moran Test if and only if the cfunc() is "1".
moran <- ifelse(moran & cfunc == "1", TRUE, FALSE)
there_are_ties <- as.logical(sum(duplicated(x)))
if(there_are_ties & ! silent) message("ties detected")
if(there_are_ties & ! ties.method == "density") {
untied <- TRUE; ITS <- 0; ix <- 1:n
while( n != length(unique(x)) ) {
ITS <- ITS + 1
for(i in 1:n) { # loop the whole data sequence, avoid R idioms here!
six <- ix[x == x[i]]; ntie <- length(six)
if(ntie == 1) next # early out and move to next point
if(ties.method == "bernstein") {
half <- ntie/2 # initial estimate of the middle
if(half != as.integer(half)) { # is odd, preserve one original value
middle.idx <- six[1] + as.integer(half)
x.in.middle <- x[middle.idx]
x[six] <- dat2bernqua(six/(n+1), x, ...)
x[middle.idx] <- x.in.middle
next;
} else { # not able to preserve an original value because of symmetry
x[six] <- dat2bernqua(six/(n+1), x, ...)
next;
}
stop("shall not be here in logic!")
}
val <- x[six[1]]; tsub <- NA
if(delta) {
del <- delta
ylo <- val - del; yhi <- val + del
tsub <- seq(ylo, yhi, by=2*del/(ntie-1))
} else {
# This ifelse() is protection against zeros
del <- ifelse(val == 0, -log10offset,
as.integer(log10(abs(val))) -log10offset)
ylo <- val - 10^del; yhi <- val + 10^del
tsub <- seq(ylo, yhi, by=2*10^del/(ntie-1))
}
if(length(tsub) != length(six)) {
stop("fatal error in algorithmic thought")
}
x[six] <- tsub
}
if(ITS > 100) {
warning("emergency break on trying to process ties, ",
"maximum iterations of 100 reached")
break
}
}
}
x <- sort(x) # final insurance on the sort
if(is.null(init.para)) {
lmr <- lmoms(x)
if(! are.lmom.valid(lmr)) {
warning("L-moments of x are not valid for initial parameters, ",
"try manual initial parameters")
return(NULL)
}
init.para <- lmom2par(lmr, type=type, ...)
if(is.null(init.para)) {
warning("could not estimate initial parameters via L-moments")
return(NULL)
}
} else if(! is.list(init.para) & is.vector(init.para)) {
init.para <- vec2par(init.para, type=type)
if(is.null(init.para)) {
warning("initial parameters given by vector are not valid for initial parameters, ",
"try other initial parameters")
return(NULL)
}
}
if(check.support) {
# plotting position a must be [0,0.5], the 1 triggers max likelihood
# and the 2 to trigger a last pass (testing the MLE) and then exit
for(a.pp.coe in c(seq(0,0.5, by=0.1), 1, 2) ) {
support <- supdist(init.para)$support
if(min(x) < support[1]) {
if(! silent) message("minimum x is < than the support of ",
"initial parameters, try alternative initial parameters")
if(a.pp.coe <= 0.5) {
if(! silent) message(" trying pwm.pp --> lmom --> a=",a.pp.coe)
lmr <- pwm2lmom(pwm.pp(x, a=a.pp.coe))
if(! are.lmom.valid(lmr)) next
init.para <- lmom2par(lmr, type=type, ...)
next
} else if(a.pp.coe == 1 & mle2par) {
if(! silent) message(" trying MLE instead for initial parameters")
init.para <- mle2par(x, type=type, ...)
if(is.null(init.para)) {
warning("attempted final spinup by subordinated call to mle2par()")
return(NULL)
}
next
} else {
warning(" giving up on automatic starting parameters")
return(NULL)
}
}
if(max(x) > support[2]) {
if(! silent) message("maximum x is > than the support of ",
"initial parameters, try alternative initial parameters")
if(a.pp.coe <= 0.5) {
if(! silent) message(" trying pwm.pp --> lmom --> para for a=",a.pp.coe)
lmr <- pwm2lmom(pwm.pp(x, a=a.pp.coe))
if(! are.lmom.valid(lmr)) next
init.para <- lmom2par(lmr, type=type, ...)
next
} else if(a.pp.coe == 1 & mle2par) {
if(! silent) message(" trying MLE instead for initial parameters")
init.para <- mle2par(x, type=type, ...)
if(is.null(init.para)) {
warning("attempted final spinup by subordinated call to mle2par()")
return(NULL)
}
next
} else {
warning(" giving up on automatic starting parameters")
return(NULL)
}
}
}
}
if(is.null(init.para)) { # finally insurance policy for error trapping
warning(" initial parameters are NULL")
return(NULL)
}
if(init.para$type != type) {
warning("distribution requested to fit does not match the type of the ",
"initial parameters")
return(NULL)
}
if(length(init.para$para) == 1) {
warning("function is not yet built for single parameter optimization")
return(NULL)
}
"afunc" <- function(para, x=NULL, n=NA, ...) {
lmomco.para <- vec2par(pretransf(para), type=type, paracheck=TRUE)
if(is.null(lmomco.para)) return(Inf) # trap if bad parameters
uu <- c(0,plmomco(x, lmomco.para),1) # padding the edges
dd <- diff(uu) # the deltas, length shrinks by one
dd[is.nan(dd)] <- 0 # this is a trap in case NaN leak through on performance of the CDF
if(ties.method == "density") {
dd[dd == 0] <- dlmomco(x[dd == 0], lmomco.para)
}
M <- sum(the.cfunc(dd))
if(! silent) message(" M=",M) # note M=0 will be shown but overrided
if(M == 0) return(Inf) # this is key to keeping the simplex alive!
return(M)
}
# print(ptransf(init.para$para))
# raw.afunc.call <- afunc(ptransf(init.para$para), x=x, n=n)
# print("RAW afunc() CALL WITH INITIAL PARAMETERS")
# print(raw.afunc.call)
# Note for some reason there is an argument name clash(?) that WHA can not get to
# the bottom of for a case where ... in the main call has say p=3 to trigger the
# generalized gamma if type="gam" so ... is not wrapped into the optim() call as
# WHA would instinctively do.
rt <- NULL
try(rt <- optim(par=ptransf(init.para$para), fn=afunc, x=x, n=n, ...), silent=silent)
if(is.null(rt)) {
warning("optim() attempt is NULL")
return(NULL)
}
if(null.on.not.converge & rt$convergence != 0) {
warning("optim() reports convergence error")
return(NULL)
}
lmomco.para <- vec2par(pretransf(rt$par), type=type)
lmomco.para$source <- "mps2par" # override vec2par()
lmomco.para$init.para <- init.para # preserve the initial parameters
M <- rt$value
if(! silent) lmomco.para$optim <- rt
moranTest <- function(M,n,p) {
moran_mean <- function(n) {
# euler <- print(-digamma(1), digits=15)
euler <- 0.5772156649015323 # Euler's constant
(n+1)*(log(n+1)+euler) - (1/2) - 1/(12*(n+1))
}
moran_variance <- function(n, minusone=FALSE) {
(n+1)*(pi^2/6 - 1) - (1/2) - 1/(6*(n+1))
}
muM <- moran_mean(n); varM <- moran_variance(n)
C1 <- muM - sqrt(varM*n/2); C2 <- sqrt(varM/(2*n))
Tstat <- (M - C1 + (p/2))/C2; Io <- (n+1)*log(n+1)
if(M < Io) {
# The check.support should keep this from ever being triggered,
# the author thinks!
warning(" Moran < Io, ",
"fitted distribution might have support that is ",
"narrower than the data, failure results in Chi-Squared")
}
ests <- data.frame(muM=muM, varM=varM, C1=C1, C2=C2, sample.size=n)
row.names(ests) <- "" # makes a print of output less confusing
f <- pchisq(Tstat, n, lower.tail=TRUE); p.value <- 1-f
zz <- round(c(Io, M, Tstat, p.value), digits=4)
names(zz) <- c( "Io", "Moran", "T(theta)", "p-value")
return(list(diagnostics=ests, statistics=zz))
}
if(there_are_ties & get.untied & ! ties.method=="density") {
lmomco.para$ties <- list(method=ties.method,
untied.pseudo.data=x, iterations=ITS)
}
if(moran) lmomco.para$MoranTest <- moranTest(M, n, length(lmomco.para$para))
return(lmomco.para)
}
# Other implementations
# fBasics --> dist-gldFit.R --> .gldFit.mps -->
# f = try(-typeFun(log(DH[DH > 0])), silent = TRUE)
# gld --> fit_fkml.R --> fit_fkml.c --> method.id == 2:
#If F[i]-F[i-1] = 0, replace by f[i-1] (ie the density at smaller observation)
# If the cfunc types are being explored:
#mean(replicate(1000, mps2par(rlmomco(100,list(type="nor", para=c(0.1,1))), type="nor", cfunc="3")$para[2]))
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.