# R/positive_mle.R In Rfast: A Collection of Efficient and Extremely Fast R Functions

#### Documented in betaprime.mlechisq.mleexp2.mlefoldnorm.mlegammamlehalfnorm.mleinvgauss.mlelindley.mlelogcauchy.mleloglogistic.mlelomax.mlemaxboltz.mlenormlog.mlepareto.mlerayleigh.mletobit.mleweibull.mle

```#[export]
betaprime.mle <- function(x, tol = 1e-09) {

n <- length(x)
slx <- sum( log(x) )
slx2 <- sum( log1p(x) )
m <- sum(x) / n    ;   s <- sum(x^2) - m^2
b <- ( m^2 + m ) / s + 2
a <- abs( m * b - m )
lik1 <- (a - 1) * slx - (a + b) * slx2 - n * lbeta(a, b)

dera <- n * digamma(a + b) - n * digamma(a) + slx - slx2
derb <- n * digamma(a + b) - n * digamma(b) - slx2
derab <- n * trigamma(a + b)
dera2 <- derab - n * trigamma(a)
derb2 <- derab - n * trigamma(b)
anew <- c(a, b) - c( derb2 * dera - derab * derb, - derab * dera + dera2 * derb ) / ( dera2 * derb2 - derab^2 )
a <- anew[1]      ;      b <- anew[2]
lik2 <- (a - 1) * slx - (a + b) * slx2 - n * lbeta(a, b)

i <- 2
while ( lik2 - lik1 > tol ) {
i <- i + 1
lik1 <- lik2
a <- anew[1]      ;      b <- anew[2]
dera <- n * digamma(a + b) - n * digamma(a) + slx - slx2
derb <- n * digamma(a + b) - n * digamma(b) - slx2
derab <- n * trigamma(a + b)
dera2 <- derab - n * trigamma(a)
derb2 <- derab - n * trigamma(b)
anew <- anew - c( derb2 * dera - derab * derb, - derab * dera + dera2 * derb ) / ( dera2 * derb2 - derab^2 )
a <- anew[1]      ;      b <- anew[2]
lik2 <- (a - 1) * slx - (a + b) * slx2 - n * lbeta(a, b)
}

names(anew) <- c("alpha", "beta")
loglik <- (a - 1) * slx - (a + b) * slx2 - n * lbeta(a, b)
list(iters = i, loglik = loglik, param = anew)
}

#[export]
chisq.mle <- function(x, tol = 1e-09) {

n <- length(x)
f1 <- 0.5 * n
f2 <- 0.5 * f1
com <-  - f1 * log(2)
slx2 <- 0.5 * sum( log(x) )
sx <- sum(x)
v <- sx / n
lik1 <- v * com - n * lgamma(0.5 * v) + 2 * (0.5 * v - 1) * slx2
der <-  com - f1 * digamma(0.5 * v) + slx2
der2 <-  - f2 * trigamma(0.5 * v)
v <- v - der / der2
i <- 2
lik2 <-  v * com - n * lgamma(0.5 * v) + 2 * (0.5 * v - 1) * slx2
while ( lik2 - lik1 > tol ) {
i <- i + 1
lik1 <- lik2
der <-  com - f1 * digamma(0.5 * v) + slx2
der2 <-  - f2 * trigamma(0.5 * v)
v <- v - der / der2
lik2 <-  v * com - n * lgamma(0.5 * v) + 2 * (0.5 * v - 1) * slx2
}

list(iters = i , loglik = lik2 - 0.5 * sx, nu = v)
}

#[export]
exp2.mle <- function(x) {
a <- min(x)
n <- length(x)
b <- sum(x)/n - a
param <- c(a, b)
names(param) <- c("a", "b")
loglik <-  - n * log(b) - 1/n
list(loglik = loglik, param = param)
}

#[export]
expmle <- function (x) {
n <- length(x)
lambda <- sum(x)/n
loglik <-  - n * log(lambda) - n
list(loglik = loglik, lambda = lambda)
}

#[export]
foldnorm.mle <- function(x, tol = 1e-09) {
n <- length(x)
m <- sum(x) / n
x2 <- x^2
sx2 <- sum(x2)
a <- sx2 / n - m^2
y <- m * x / a
tanhy <- tanh(y)
coshy2 <- 1 / cosh(y)^2
derm <-  - n * m / a + sum(x * tanhy ) / a
ders <-  - n / 2 / a + (sx2 + n * m^2) / 2 / a^2 - sum( y * tanhy ) / a
derm2 <-  - n/ a + sum(x2 * coshy2 ) / a^2
ders2 <- n / 2 / a^2 - (sx2 + n * m^2) / a^3  + 2 * sum( y * tanhy ) / a^2 + sum(y^2 * coshy2 ) / a^2
derms <- n * m / a^2 - sum(x * tanhy ) / a^2 - sum(y * x * coshy2 ) / a^2
aold <- c(m, a)
anew <- aold - c( ders2 * derm - derms * ders, - derms * derm + derm2 * ders) / (derm2 * ders2 - derms^2)

i <- 2
while ( sum( abs(anew - aold) ) > tol ) {
i <- i + 1
m <- anew[1]    ;   a <- anew[2]
aold <- anew
y <- m * x / a
tanhy <- tanh(y)
coshy2 <- 1 / cosh(y)^2
derm <-  - n * m / a + sum(x * tanhy ) / a
ders <-  - n / 2 / a + (sx2 + n * m^2) / 2 / a^2 - sum( y * tanhy ) / a
derm2 <-  - n/ a + sum(x2 * coshy2 ) / a^2
ders2 <- n / 2 / a^2 - (sx2 + n * m^2) / a^3  + 2 * sum( y * tanhy ) / a^2 + sum(y^2 * coshy2 ) / a^2
derms <- n * m / a^2 - sum(x * tanhy ) / a^2 - sum(y * x * coshy2 ) / a^2
anew <- aold - c( ders2 * derm - derms * ders, - derms * derm + derm2 * ders) / (derm2 * ders2 - derms^2)
}

m <- anew[1]   ;  a <- anew[2]
names(anew) <- c("mean", "sigma squared")
loglik <- n/2 * log( 2 / pi / a) - 0.5 * n * m^2 / a - 0.5 * sum(x2) / a + sum( log( cosh( y ) ) )
list(iters = i, loglik = loglik, param = anew)
}

#[export]
gammamle <- function(x, tol = 1e-09) {
n <- length(x)
m <- sum(x)/n
slx <- sum( log(x) ) / n
s <- log(m) - slx
a1 <- 3 - s + sqrt( (s-3)^2 + 24 * s )
a1 <- a1 / (12 * s)
a2 <- a1 - ( log(a1) - digamma(a1) - s) / (1/a1 - trigamma(a1) )
i <- 2
while ( abs(a2 - a1) > tol) {
i <- i + 1
a1 <- a2
a2 <- a1 - ( log(a1) - digamma(a1) - s) / (1/a1 - trigamma(a1) )
}
b <- a2 / m
loglik <-  - b * n * m + (a2 - 1) * n * slx + n * a2 * log(b) - n *
lgamma(a2)
param <- c(a2, b)
names(param) <- c("shape", "rate")
list(iters = i, loglik = loglik, param = param)
}

#old_gammamle <- function(x, tol = 1e-09) {
#  n <- length(x)
#  sx <- sum(x)
#  sx2 <- sum(x^2)
#  slx <- sum( log(x) )
# b <-  sx / sx2     ;     a <- b * sx / n
# dera <- slx + n * log(b) - n * digamma(a)
# dera2 <-  - n * trigamma(a)
# derb2 <-  - n * a / b^2
# derb <-  - sx  - b * derb2
# derab <-  n / b
# aold <- c(a, b)
# anew <- aold - c( derb2 * dera - derab * derb, - derab * dera + dera2 * derb ) / ( dera2 * derb2 - derab^2 )
# i <- 2
# while ( sum( abs(anew - aold) ) > tol ) {
# i <- i + 1
# aold <- anew
# a <- anew[1]     ;      b <- anew[2]
# dera <- slx + n * log(b) - n * digamma(a)
# dera2 <-  - n * trigamma(a)
# derb2 <-  - n * a / b^2
# derb <-  - sx  - b * derb2
# derab <-  n / b
# anew <- aold - c( derb2 * dera - derab * derb, - derab * dera + dera2 * derb ) / ( dera2 * derb2 - derab^2 )
# }
# a <- anew[1]    ;     b <- anew[2]
# loglik <-  - b * sx + (a - 1) * slx + n * a * log(b) - n * lgamma(a)
# names(anew) <- c("shape", "scale")
# list(iters = i, loglik = loglik, param = anew)
# }

#[export]
halfnorm.mle <- function(x) {
n <- length(x)
s <- sqrt( sum(x^2) / n )
loglik <- n / 2 * log(2 / s / pi) - n / 2
list(loglik = loglik, sigma.squared = s)
}

#[export]
invgauss.mle <- function(x) {
n <- length(x)
sx <- sum(x)
sx2 <- sum(1/x)
m <- sx / n
lambda <- 1 / ( sx2 / n - 1 / m )
loglik <- n * 0.5 * log(lambda / 2 / pi) - 1.5 * sum( log(x) ) - lambda / 2 / m^2 * (-sx + m^2 * sx2)
param <- c(m, lambda)
names(param) <- c("mu", "lambda")
list(loglik = loglik, param = param)
}

#[export]
lindley.mle <- function(x) {
n <- length(x)
sx <- sum(x)
a <- sx/n
b <- a - 1
delta <- b^2 + 8 * a
theta <-  0.5 * ( - b + sqrt(delta) ) / a
loglik <- 2 * n * log(theta) - n * log(1 + theta) + sum( log1p(x) ) - theta * sx
list(loglik = loglik, theta = theta)
}

#[export]
logcauchy.mle <- function(x, tol = 1e-09) {
x <- log(x)
a <- Rfast::cauchy.mle(x, tol = tol)\$param
m <- a[1]  ;  s <- a[2]
loglik <- length(x) * log(s/pi) - sum(x) - sum( log( (x - m)^2 + s^2) )
list(loglik = loglik, param = a)
}

#[export]
loglogistic.mle <- function(x, tol = 1e-09) {
y <- log(x)
n <- length(x)
param <- Rfast::logistic.mle(y, tol = tol)\$param
a <- exp(param[1])   ;  b <- 1/param[2]
loglik <- n * log(b/a) + (b - 1) * sum(y) - n * (b - 1) * log(a) - 2 * sum( log1p( (x/a)^b ) )
list(loglik = loglik, param = c(a, b) )
}

#[export]
lognorm.mle <- function (x) {
n <- length(x)
x <- log(x)
sx <- sum(x)
m <- sx/n
s <- sum(x^2)/n - m^2
loglik <-  - 0.5 * n * ( log(2 * pi * s) + 1 ) - sx
param <- c(m, s)
names(param) <- c("mean", "variance")
list(loglik = loglik, param = param)
}

#[export]
lomax.mle <- function(x, tol = 1e-09) {
n <- length(x)
x2 <- x^2
m <- sum(x) / n
s2 <- sum(x2) / n - m^2
expa <- abs( 2 * s2 / (s2 - m^2) )
explam <- abs( expa - 1 ) * m
dera2 <-  - expa * sum( log1p( x / explam) )
dera <-  n + dera2
com <- sum( x / (explam + x) )
derlama <- expa * com
derlam <-  - n + derlama + com
derlam2 <-  - ( expa + 1)* explam * sum( x / (explam + x)^2 )
aold <- log( c(expa, explam) )
anew <- aold - c( derlam2 * dera - derlama * derlam, - derlama * dera + dera2 * derlam ) / ( dera2 * derlam2 - derlama^2 )

i <- 2
while ( sum( abs(anew - aold) ) > tol ) {
i <- i + 1
aold <- anew
expa <- exp(anew[1])     ;      explam <- exp(anew[2])
dera2 <-  - expa * sum( log1p( x / explam) )
dera <-  n + dera2
com <- sum( x / (explam + x) )
derlama <- expa * com
derlam <-  - n + derlama + com
derlam2 <-  - ( expa + 1) * explam * sum( x / (explam + x)^2 )
anew <- aold - c( derlam2 * dera - derlama * derlam, - derlama * dera + dera2 * derlam ) / ( dera2 * derlam2 - derlama^2 )
}

a <- exp(anew[1])    ;     lam <- exp(anew[2])
loglik <- n * log(a / lam) - (a + 1) * sum( log1p(x / lam) )
names(anew) <- c("shape", "scale")
list( iters = i, loglik = loglik, param = c(a, lam) )
}

#[export]
maxboltz.mle <- function(x) {
n <- length(x)
a <- sqrt( sum(x^2) / 3 / n )
loglik <- n/2 * log(2 / pi) +  2 * sum( log(x) ) - 1.5 * n - 3 * n * log(a)
list(loglik = loglik, a = a)
}

#[export]
normlog.mle <- function(x) {
n <- length(x)
mx <- sum(x) /n
m <- log( mx )
sigma <- sum(x^2)/n - mx^2
loglik <-  - 0.5 * n * log(2 * pi * sigma) - 0.5 * n
param <- c(mx, m, sigma, sigma * n / (n - 1) )
names(param) <- c("exp_mu", "mu", "biased variance", "unbiased variance")
list(loglik = loglik, param = param)
}

#[export]
pareto.mle <- function(x) {
n <- length(x)
xm <- min(x)
com <- n * log(xm)
slx <- sum( log(x) )
a <- n / ( slx - com )
param <- c(xm, a)
names(param) <- c("scale", "shape")
loglik <- n * log(a) + a * com - (a + 1) * slx
list(loglik = loglik, param = param)
}

#[export]
rayleigh.mle <- function(x) {
n <- length(x)
sigma <- 0.5 * sum(x^2) / n
loglik <- sum( log(x) ) - n * log(sigma) - n
list(loglik = loglik, sigma = sigma)
}

#[export]
tobit.mle <- function(y, tol = 1e-09) {
y1 <- y[y >0]  ;  n1 <- length(y1)
n <- length(y)
n0 <- n - n1
sy12 <- sum(y1^2)
m <- mean(y)   ;  s <- sqrt( sy12/n - m^2 )
sy1 <- n * m
com <- dnorm(m, 0, s) / pnorm(-m/s)
derm <- (sy1 - n1 * m)/s^2 - n0 * com
derm2 <-  -n1/s^2 - n0 * ( -m /s^2 * com + com^2 )
ders <-  - n1 + (sy12 - 2 * m * sy1 + n1 * m^2)/s^2 + n0 * m * com
ders2 <-  - 2 * (sy12 - 2 * m * sy1 + n1 * m^2)/s^2 + n0 * m * ( - com + m^2/s^2 * com - com^2 * m )
derms <-  - 2 * (sy1 - n1 * m)/s^2 - n0 * ( - com + m^2/s^2 * com - com^2 * m )
aold <- c(m, log(s))
anew <- aold - c( ders2 * derm - derms * ders, - derms * derm + derm2 * ders ) / ( derm2 * ders2 - derms^2 )
i <- 2
while ( sum( abs(aold - anew) ) > tol ) {
i <- i + 1
aold <- anew
m <- anew[1]     ;    s <- exp( anew[2] )
com <- dnorm(m, 0, s) / pnorm(-m/s, 0, 1)
derm <- (sy1 - n1 * m)/s^2 - n0 * com
derm2 <-  - n1/s^2 - n0 * ( -m /s^2 * com + com^2 )
ders <-  - n1 + (sy12 - 2 * m * sy1 + n1 * m^2)/s^2 + n0 * m * com
ders2 <-  - 2 * (sy12 - 2 * m * sy1 + n1 * m^2)/s^2 + n0 * m * ( - com + m^2/s^2 * com - com^2 * m )
derms <-  - 2 * (sy1 - n1 * m)/s^2 + n0 * ( com - m^2/s^2 * com + com^2 * m )
anew <- aold - c( ders2 * derm - derms * ders, - derms * derm + derm2 * ders ) / ( derm2 * ders2 - derms^2 )
}
s <- exp(anew[2])
loglik <-  - 0.5 * n1 * log(2 * pi * s^2) - 0.5 * ( sy12 - 2 * m * sy1 + n1 * m^2 ) / s^2 + n0 * pnorm(-m/s, log.p = TRUE)
param <- c(anew[1], s)
names(param) <- c("location", "scale")
list(iters = i, loglik = loglik, param = param)
}

#[export]
weibull.mle <- function(x, tol = 1e-09, maxiters = 100) {
mod <- .Call(Rfast_weibull_mle, x, tol, maxiters)
names(mod\$param) <- c("shape", "scale")
mod
}
```

## Try the Rfast package in your browser

Any scripts or data that you put into this service are public.

Rfast documentation built on Feb. 16, 2023, 5:09 p.m.