    ## 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
    }