## Probability density function for the Generalised Normal Laplace distribution
dgnl <- function(x, mu = 0, sigma = 1, alpha = 1, beta = 1, rho = 1,
param = c(mu, sigma, alpha, beta, rho)) {
## check parameters
parResult <- gnlCheckPars(param)
case <- parResult$case
errMessage <- parResult$errMessage
if (case == "error")
stop(errMessage)
mu <- param[1]
sigma <- param[2]
alpha <- param[3]
beta <- param[4]
rho <- param[5]
## Shifting by mu
x <- x - mu
## Initialising result vector
pdfValues <- rep(0, length(x))
## Because 'integrate' doesn't take vectors as input, we need to iterate over
## x to evaluate densities
for (i in 1:length(x)) {
## Modified characteristic function. Includes minor calculation regarding
## complex numbers to ensure the function returns a real number
chfn <- function(s) {
result <- (alpha * beta * exp(-((sigma^2 * s^2) / 2))) /
(complex(real = alpha, imaginary = -s) *
complex(real = beta, imaginary = s))
result <- result^rho ## Scaling result by rho
r <- Mod(result)
theta <- Arg(result)
r * cos(theta - (s * x[i]))
}
## Integrating modified characteristic function
pdfValues[i] <- (1 / pi) * integrate(chfn, 0, Inf)$value
}
## Returning vector of densities
pdfValues
}