# Copyright (c) 2015, Roberto Garrappa
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
# * Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
# * Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in
# the documentation and/or other materials provided with the distribution
# * Neither the name of the Department of Mathematics - University of Bari -
# Italy nor the names of its contributors may be used to endorse or promote
# products derived from this software without specific prior written
# permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
# The code below is adapted from:
# https://au.mathworks.com/matlabcentral/fileexchange/48154-the-mittag-leffler-function
#' @name MittagLeffleR
NULL
#' @rdname MittagLeffleR
#' @export
#' @param z The argument (real-valued)
#' @param a,b,g Parameters of the Mittag-Leffler distribution; see Garrappa
#' @return \code{mlf} returns the value of the Mittag-Leffler function.
#' @examples
#' mlf(2,0.7)
mlf <- function(z,a,b=1,g=1) {
#check number of inputs and handle errors
if (nargs() < 2) {
#handle error
}
#check inputs values and handle errors
#target precision
log_epsilon <- log(10^(-15))
#LT inversion for each element of z
E <- numeric(length(z))
for (i in 1:length(z)) {
if (abs(z[i]) < 10^(-15) ) {
E[i] <- 1/gamma(b)
}
else {
E[i] <- LTinversion(1,z[i],a,b,g,log_epsilon)
}
}
return(E)
}
LTinversion <- function(t, lambda, a, b, g, log_epsilon) {
theta <- Arg(lambda)
kmin <- ceiling(-a/2 - theta/2/pi)
kmax <- floor(a/2 - theta/2/pi)
if (kmin > kmax) {
k_vett <- c()
}
else {
k_vett <- kmin:kmax
}
s_star <- abs(lambda)^(1/a) * exp(1i*(theta+2*k_vett*pi)/a)
phi_s_star <- (Re(s_star)+abs(s_star))/2
phi_s_star <- sort(phi_s_star, index.return=TRUE)$x
index_s_star <- sort(phi_s_star, index.return=TRUE)$ix
s_star <- s_star[index_s_star]
index_save <- phi_s_star > 10^(-15)
s_star <- s_star[index_save]
phi_s_star <- phi_s_star[index_save]
s_star <- c(0, s_star)
phi_s_star <- c(0, phi_s_star)
J1 <- length(s_star)
J <- J1 - 1
p <- c( max(0,-2*(a*g-b+1)) , rep(1,J)*g )
q <- c( rep(1,J)*g , +Inf)
phi_s_star <- c(phi_s_star, +Inf)
admissible_regions <- which( (phi_s_star[1:length(phi_s_star)-1] < (log_epsilon - log(.Machine$double.eps))/t) & (phi_s_star[1:length(phi_s_star)-1] < phi_s_star[2:length(phi_s_star)]))
JJ1 <- admissible_regions[length(admissible_regions)]
mu_vett <- rep(1,JJ1)*Inf
N_vett <- rep(1,JJ1)*Inf
h_vett <- rep(1,JJ1)*Inf
find_region <- 0
while (!find_region) {
for (j1 in admissible_regions) {
if (j1 < J1) {
temp2 <- OptimalParam_RB(t,phi_s_star[j1],phi_s_star[j1+1],p[j1],q[j1],log_epsilon)
muj <- temp2[1]
hj <- temp2[2]
Nj <- temp2[3]
}
else {
temp3 <- OptimalParam_RU(t,phi_s_star[j1],p[j1],log_epsilon)
muj <- temp3[1]
hj <- temp3[2]
Nj <- temp3[3]
}
mu_vett[j1] <- muj
h_vett[j1] <- hj
N_vett[j1] <- Nj
}
if (min(N_vett) > 200) {
log_epsilon <- log_epsilon +log(10)
}
else {
find_region <- 1
}
}
N <- min(N_vett)
iN <- which.min(N_vett)
mu <- mu_vett[iN]
h <- h_vett[iN]
k <- -N:N
u <- h*k
z <- mu*(1i*u+1)^2
zd <- -2*mu*u + 2*mu*1i
zexp <- exp(z*t)
F <- z^(a*g-b)/(z^a - lambda)^g*zd
S <- zexp*F
Integral <- h*sum(S)/2/pi/1i
if ((iN+1) > length(s_star)) {
ss_star <- c()
}
else {
ss_star <- s_star[(iN+1):length(s_star)]
}
Residues <- sum(1/a*(ss_star)^(1-b)*exp(t*ss_star))
E <- Integral + Residues
if (!is.complex(lambda)) {
E <- Re(E)
}
return(E)
}
OptimalParam_RB <- function(t, phi_s_star_j, phi_s_star_j1, pj, qj, log_epsilon) {
log_eps <- -36.043653389117154 #log(eps)
fac <- 1.01
conservative_error_analysis <- 0
f_max <- exp(log_epsilon - log_eps)
sq_phi_star_j <- sqrt(phi_s_star_j)
threshold <- 2*sqrt((log_epsilon - log_eps)/t) ;
sq_phi_star_j1 <- min(sqrt(phi_s_star_j1), threshold - sq_phi_star_j) ;
if ((pj < 10^(-14)) && (qj < 10^(-14))) {
sq_phibar_star_j <- sq_phi_star_j
sq_phibar_star_j1 <- sq_phi_star_j1
adm_region <- 1
}
if ((pj < 10^(-14)) && (qj >= 10^(-14))) {
sq_phibar_star_j <- sq_phi_star_j
if (sq_phi_star_j > 0) {
f_min <- fac*(sq_phi_star_j/(sq_phi_star_j1-sq_phi_star_j))^qj
}
else {
f_min <- fac
}
if (f_min < f_max) {
f_bar <- f_min + f_min/f_max*(f_max-f_min)
fq <- f_bar^(-1/qj)
sq_phibar_star_j1 <- (2*sq_phi_star_j1 - fq*sq_phi_star_j)/(2+fq)
adm_region <- 1
}
else {
adm_region <- 0
}
}
if ((pj >= 10^(-14)) && (qj < 10^(-14))) {
sq_phibar_star_j1 <- sq_phi_star_j1
f_min = fac*(sq_phi_star_j1/(sq_phi_star_j1 - sq_phi_star_j))^pj
if (f_min < f_max) {
f_bar <- f_min + f_min/f_max*(f_max - f_min)
fp <- f_bar^(-1/pj)
sq_phibar_star_j <- (2*sq_phi_star_j+fp*sq_phi_star_j1)/(2-fp)
adm_region <- 1
}
else {
adm_region <- 0
}
}
if ((pj >= 10^(-14)) && (qj >= 10^(-14))) {
f_min = fac*(sq_phi_star_j+sq_phi_star_j1)/(sq_phi_star_j1-sq_phi_star_j)^max(pj,qj)
if (f_min < f_max) {
f_min <- max(f_min,1.5)
f_bar <- f_min + f_min/f_max*(f_max-f_min)
fp <- f_bar^(-1/pj)
fq <- f_bar^(-1/qj)
if (!conservative_error_analysis) {
w <- -phi_s_star_j1*t/log_epsilon
}
else {
w <- -2*phi_s_star_j1*t/(log_epsilon-phi_s_star_j1*t)
}
den <- 2+w - (1+w)*fp + fq
sq_phibar_star_j <- ((2+w+fq)*sq_phi_star_j + fp*sq_phi_star_j1)/den
sq_phibar_star_j1 <- (-(1+w)*fq*sq_phi_star_j + (2+w-(1+w)*fp)*sq_phi_star_j1)/den
adm_region <- 1
}
else {
adm_region <- 0
}
}
if (adm_region) {
log_epsilon <- log_epsilon - log(f_bar)
if (!conservative_error_analysis) {
w <- -sq_phibar_star_j1^2*t/log_epsilon
}
else {
w = -2*sq_phibar_star_j1^2*t/(log_epsilon-sq_phibar_star_j1^2*t)
}
muj <- (((1+w)*sq_phibar_star_j + sq_phibar_star_j1)/(2+w))^2
hj <- -2*pi/log_epsilon*(sq_phibar_star_j1-sq_phibar_star_j)/((1+w)*sq_phibar_star_j + sq_phibar_star_j1)
Nj <- ceiling(sqrt(1-log_epsilon/t/muj)/hj)
}
else {
muj <- 0
hj <- 0
Nj <- +Inf
}
out <- c(muj, hj, Nj)
return(out)
}
OptimalParam_RU <- function(t, phi_s_star_j, pj, log_epsilon) {
sq_phi_s_star_j <- sqrt(phi_s_star_j)
if (phi_s_star_j > 0) {
phibar_star_j <- phi_s_star_j*1.01
}
else {
phibar_star_j <- 0.01
}
sq_phibar_star_j <- sqrt(phibar_star_j)
f_min <- 1
f_max <- 10
f_tar <- 5
stopp <- 0
while (!stopp) {
phi_t <- phibar_star_j*t
log_eps_phi_t <- log_epsilon/phi_t
Nj <- ceiling(phi_t/pi*(1 - 3*log_eps_phi_t/2 + sqrt(1-2*log_eps_phi_t)))
A <- pi*Nj/phi_t
sq_muj <- sq_phibar_star_j*abs(4-A)/abs(7-sqrt(1+12*A))
fbar <- ((sq_phibar_star_j-sq_phi_s_star_j)/sq_muj)^(-pj)
stopp <- (pj < 10^(-14)) || (f_min < fbar && fbar < f_max)
if (!stopp) {
sq_phibar_star_j <- f_tar^(-1/pj)*sq_muj + sq_phi_s_star_j
phibar_star_j <- sq_phibar_star_j^2
}
}
muj <- sq_muj^2
hj <- (-3*A - 2 + 2*sqrt(1+12*A))/(4-A)/Nj
log_eps <- log(.Machine$double.eps)
threshold <- (log_epsilon - log_eps)/t
if (muj > threshold) {
if (abs(pj) < 10^(-14)) {
Q <- 0
}
else {
Q <- f_tar^(-1/pj)*sqrt(muj)
}
phibar_star_j <- (Q + sqrt(phi_s_star_j))^2
if (phibar_star_j < threshold) {
w <- sqrt(log_eps/(log_eps-log_epsilon))
u <- sqrt(-phibar_star_j*t/log_eps)
muj <- threshold
Nj <- ceiling(w*log_epsilon/2/pi/(u*w-1))
hj <- sqrt(log_eps/(log_eps - log_epsilon))/Nj
}
else {
Nj <- +Inf
hj = 0
}
}
out <- c(muj, hj, Nj)
return(out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.