R/support.R In rmetalog: The Metalog Distribution

Defines functions pdfMetalog_densitynewtons_method_metalogdiffMatMetalogquantileMetalogpdfMetalogMLprobs

# Supporting functions called inside the metalog function call

# Build the quantiles through a base function
MLprobs <- function(x, step_len) {
if (class(x) != 'numeric') {
return(print('Error: input must be a numeric vector!'))
}

# Need to sort the dataframe for the quantiles to work out
l <- length(x)
x <- sort(x)
x <- as.data.frame(x)

# Calculate the likelihood as an interpolation
x\$probs <- 0
for (i in 1:l) {
if (i == 1) {
x\$probs[i] <- (0.5 / l)
}
else{
x\$probs[i] <- (x\$probs[i - 1] + (1 / l))
}
}

# If the data is very long we down convert to a smaller but representative
# vector using the step_len default is 0.01 which is a 109 element vector with
# fine values in the tail (tailstep)
if (nrow(x) > 100) {
y <- seq(step_len, (1 - step_len), step_len)
tailstep <- (step_len / 10)
y <- c(seq(tailstep, (min(y) - tailstep), tailstep),
y,
seq((max(y) + tailstep), (max(y) + tailstep * 9), tailstep))
x_new <- stats::quantile(x[, 1], probs = y)
x <- as.data.frame(x_new)
x\$probs <- y
}

return(x)
}

pdfMetalog <- function(a,
y,
t,
bounds = c(),
boundedness = 'u') {
d <- y * (1 - y)
f <- (y - 0.5)
l <- log(y / (1 - y))

# Initiate pdf

# For the first three terms
x <- (a / d)
if (a != 0) {
x <- x + a * ((f / d) + l)
}

# For the fourth term
if (t > 3) {
x <- x + a
}

# Initalize some counting variables
e <- 1
o <- 1

# For all other terms greater than 4
if (t > 4) {
for (i in 5:t) {
if (i %% 2 != 0) {
# iff odd
x <- x + ((o + 1) * a[i] * f^o)
o <- o + 1
}

if (i %% 2 == 0) {
# iff even
x <- x + a[i] * (((f^(e + 1)) / d) + (e + 1) * (f^e) * l)
e <- e + 1
}
}
}

# Some change of variables here for boundedness

x <- (x^(-1))

if (boundedness != 'u') {
M <- quantileMetalog(a, y, t, bounds = bounds, boundedness = 'u')
}

if (boundedness == 'sl') {
x <- x * exp(-M)
}

if (boundedness == 'su') {
x <- x * exp(M)
}

if (boundedness == 'b') {
x <- (x * (1 + exp(M))^2) / ((bounds - bounds) * exp(M))
}

return(x)
}

# Quantile function
quantileMetalog <- function(a,
y,
t,
bounds = c(),
boundedness = 'u') {
# Some values for calculation
f <- (y - 0.5)
l <- log(y / (1 - y))

# For the first three terms
x <- a + a * l + a * f * l

# For the fourth term
if (t > 3) {
x <- x + a * f
}

# Some tracking variables
o <- 2
e <- 2

# For all other terms greater than 4
if (t > 4) {
for (i in 5:t) {
if (i %% 2 == 0) {
x <- x + a[i] * f^e * l
e <- e + 1
}

if (i %% 2 != 0) {
x <- x + a[i] * f^o
o <- o + 1
}
}
}

if (boundedness == 'sl') {
x <- bounds + exp(x)
}

if (boundedness == 'su') {
x <- bounds - exp(-x)
}

if (boundedness == 'b') {
x <- (bounds + bounds * exp(x)) / (1 + exp(x))
}

return(x)
}

# Function for returning the matrix of differentiation terms
diffMatMetalog <- function(term_limit, step_len) {
y <- seq(step_len, (1 - step_len), step_len)
Diff <- data.frame()

for (i in 1:length(y)) {
d <- y[i] * (1 - y[i])
f <- (y[i] - 0.5)
l <- log(y[i] / (1 - y[i]))

# Initiate pdf
diffVector <- 0

# For the first three terms
x <- (1 / d)
diffVector <- c(diffVector, x)

if (term_limit > 2) {
diffVector <- c(diffVector, ((f / d) + l))
}

# For the fourth term
if (term_limit > 3) {
diffVector <- c(diffVector, 1)
}

# Initalize some counting variables
e <- 1
o <- 1

# For all other terms greater than 4
if (term_limit > 4) {
for (i in 5:term_limit) {
if (i %% 2 != 0) {
# iff odd
diffVector <- c(diffVector, ((o + 1) * f^o))
o <- o + 1

}
if (i %% 2 == 0) {
# iff even
diffVector <- c(diffVector,
(((f^(e + 1)) / d) + (e + 1) * (f^e) * l))
e <- e + 1
}
}
}

Diff <- rbind(Diff, diffVector)
}

Diff <- as.matrix(Diff)
Diff_neg = -(Diff)
new_Diff <- matrix(c(Diff[, 1], Diff_neg[, 1]), ncol = 2)

for (c in 2:length(Diff[1, ])) {
new_Diff <- cbind(new_Diff, Diff[, c])
new_Diff <- cbind(new_Diff, Diff_neg[, c])
}

return(new_Diff)
}

newtons_method_metalog <- function(m,q,term){
#a simple newtons method application
alpha_step<-0.01
err<-0.0000001
temp_err<-0.1
y_now<-0.5
avec<-paste0('a',term)
a<-m\$A[,`avec`]
i<-1
while(temp_err>err){
frist_function<-(quantileMetalog(a,y_now,term,m\$params\$bounds,m\$params\$boundedness)-q)
derv_function<-pdfMetalog(a,y_now,term,m\$params\$bounds,m\$params\$boundedness)
y_next<-y_now-alpha_step*(frist_function*derv_function)
temp_err<-abs((y_next-y_now))
if(y_next>1){
y_next<-0.99999
}
if(y_next<0){
y_next<-0.000001
}
y_now<-y_next
i<-i+1
if(i>10000){
stop(paste0('Approximation taking too long, quantile value: ',q,' is to far from distribution median. Try plot() to see distribution.'))
}
}

return(y_now)
}

pdfMetalog_density <- function(m,t,y) {

avec<-paste0('a',t)
a<-m\$A[,`avec`]
bounds<-m\$params\$bounds
boundedness<-m\$params\$boundedness

d <- y * (1 - y)
f <- (y - 0.5)
l <- log(y / (1 - y))

# Initiate pdf

# For the first three terms
x <- (a / d)
if (a != 0) {
x <- x + a * ((f / d) + l)
}

# For the fourth term
if (t > 3) {
x <- x + a
}

# Initalize some counting variables
e <- 1
o <- 1

# For all other terms greater than 4
if (t > 4) {
for (i in 5:t) {
if (i %% 2 != 0) {
# iff odd
x <- x + ((o + 1) * a[i] * f^o)
o <- o + 1
}

if (i %% 2 == 0) {
# iff even
x <- x + a[i] * (((f^(e + 1)) / d) + (e + 1) * (f^e) * l)
e <- e + 1
}
}
}

# Some change of variables here for boundedness

x <- (x^(-1))

if (boundedness != 'u') {
M <- quantileMetalog(a, y, t, bounds = bounds, boundedness = 'u')
}

if (boundedness == 'sl') {
x <- x * exp(-M)
}

if (boundedness == 'su') {
x <- x * exp(M)
}

if (boundedness == 'b') {
x <- (x * (1 + exp(M))^2) / ((bounds - bounds) * exp(M))
}

return(x)
}

Try the rmetalog package in your browser

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

rmetalog documentation built on Jan. 26, 2021, 1:07 a.m.