R/functions.R

Defines functions mleALD minfunc likALD absALD kurtALD skewALD varALD meanALD momentALD rALD qALD pALD dALD

Documented in absALD dALD kurtALD likALD meanALD mleALD momentALD pALD qALD rALD skewALD varALD

##########################################################################################
#ASSYMETRIC LAPLACE DISTRIBUTION (d,p,q and r functions)
##########################################################################################

dALD = function(y,mu=0,sigma=1,p=0.5)
{
  if(length(y) == 0) stop("y must be provided.")
  if(sum(y[is.na(y)==TRUE]) > 0) stop("There are some NA's values in y.")
  if(sigma <= 0) stop("sigma must be a positive number.")
  if(p >= 1 | p <= 0) stop("p must be a real number in (0,1).")
  if(abs(mu) ==Inf) stop("mu must be a finite real number.")
  
  densi = ifelse(test=y<mu,yes=(p*(1-p)/sigma)*exp((1-p)*(y-mu)/sigma),no=(p*(1-p)/sigma)*exp(-(p)*(y-mu)/sigma))
  return(densi)
}

pALD = function(q,mu=0,sigma=1,p=0.5,lower.tail=TRUE)
{
  if(length(q) == 0) stop("q must be provided.")
  if(sigma <= 0) stop("sigma must be a positive number.")
  if(p >= 1 | p <= 0) stop("p must be a real number in (0,1).")
  if(abs(mu) ==Inf) stop("mu must be a finite real number.")
  if(lower.tail != TRUE && lower.tail != FALSE) stop("lower.tail must be TRUE or FALSE.")
  
  acum = ifelse(test=q<mu,yes=p*exp((1-p)*(q-mu)/sigma),no=1-(1-p)*exp(-(p)*(q-mu)/sigma))
  ifelse(test=lower.tail == TRUE,yes=return(acum),no=return(1-acum))
}

qALD = function(prob,mu=0,sigma=1,p=0.5,lower.tail=TRUE)
{
  if(length(prob) == 0) stop("prob must be provided.")
  if(sigma <= 0) stop("sigma must be a positive number.")
  if(p >= 1 | p <= 0) stop("p must be a real number in (0,1).")
  if(abs(mu) ==Inf) stop("mu must be a finite real number.")
  if(lower.tail != TRUE && lower.tail != FALSE) stop("lower.tail must be TRUE or FALSE.")
  if(sum(prob > 1 | prob < 0) > 0) stop("All elements of prob must be real numbers in [0,1].")
  
  q = sapply(X=prob,FUN=function(prob,mu,sigma,p){ifelse(test=lower.tail == TRUE,yes=prob,no=(1-prob));
                                                  ifelse(test=prob<p,yes=mu+((sigma*log(prob/p))/(1-p)),no=mu-((sigma*log((1-prob)/(1-p)))/p))},mu=mu,sigma=sigma,p=p)
  return(q)
}

rALD = function(n,mu=0,sigma=1,p=0.5)
{
  if(length(n) == 0) stop("The sample size n must be provided.")
  if(n <= 0 || n%%1 !=0) stop("The sample size n must be a positive integer value.")
  if(sigma <= 0) stop("sigma must be a positive number.")
  if(p >= 1 | p <= 0) stop("p must be a real number in (0,1).")
  if(abs(mu) ==Inf) stop("mu must be a finite real number.")
  
  u = runif(n)
  r = sapply(X=u,FUN=qALD,mu=mu,sigma=sigma,p=p)
  return(r)
}


momentALD = function(k=1,mu=0,sigma=1,p=0.5)
{
  if(k <= 0 || k%%1 !=0) stop("The moment must be a positive integer value.")
  if(sigma <= 0) stop("sigma must be a positive number.")
  if(p >= 1 | p <= 0) stop("p must be a real number in (0,1).")
  if(abs(mu) ==Inf) stop("mu must be a finite real number.")
  
  
  moment = factorial(k)*(sigma^k)*p*(1-p)*(1/(p^(k+1)) + ((-1)^k)/((1-p)^(k+1)))
  return(moment)
}

meanALD = function(mu=0,sigma=1,p=0.5)
{
  mean = momentALD(k=1,mu=mu,sigma=sigma,p=p) + mu
  return(mean)
}

varALD = function(mu=0,sigma=1,p=0.5)
{
  var = (momentALD(k=2,mu=mu,sigma=sigma,p=p) - (momentALD(k=1,mu=mu,sigma=sigma,p=p)^2))
  return(var)
}


skewALD = function(mu=0,sigma=1,p=0.5)
{
  skew = (momentALD(k=3,mu=mu,sigma=sigma,p=p) - 
            3*meanALD(mu,sigma,p)*varALD(mu,sigma,p) - 
            meanALD(mu,sigma,p)^3)/(varALD(mu,sigma,p)^(3/2))
  return(skew)
}

kurtALD = function(mu=0,sigma=1,p=0.5)
{
  kurt = (momentALD(k=4,mu=mu,sigma=sigma,p=p) - 
            4*meanALD(mu,sigma,p)*momentALD(k=3,mu=mu,sigma=sigma,p=p) + 
            6*(meanALD(mu,sigma,p)^2)*momentALD(k=2,mu=mu,sigma=sigma,p=p) -
            3*(meanALD(mu,sigma,p)^4))/(varALD(mu,sigma,p)^2)
  return(kurt)
}

absALD = function(sigma=1,p=0.5)
{
  if(sigma <= 0) stop("sigma must be a positive number.")
  if(p >= 1 | p <= 0) stop("p must be a real number in (0,1).")
  
  absmoment = (sigma)*((1 - 2*p + 2*p^2)/(p*(1-p)))
  return(absmoment)
}

# momentALD(k=3,mu=0,sigma=1,p=0.5)
# meanALD(mu=0,sigma=1,p=0.5)
# varALD(mu=0,sigma=1,p=0.5)
# skewALD(mu=0,sigma=1,p=0.5)
# kurtALD(mu=0,sigma=1,p=0.5)
# absALD(sigma=1,p=0.5)
# 
# likALD(y,mu=0,sigma=1,p=0.5,loglik = F)


likALD = function(y,mu=0,sigma=1,p=0.5,loglik=TRUE)
{
  if(loglik != TRUE && loglik != FALSE) stop("log must be TRUE or FALSE.")
  ifelse(test=loglik == FALSE,yes=return(prod(dALD(y,mu,sigma,p=p))),no=return(sum(log(dALD(y,mu,sigma,p=p)))))
}


minfunc = function(y,mu,p)
{
  if(length(y) == 0) stop("y must be provided.")
  if(sum(y[is.na(y)==TRUE]) > 0) stop("There are some NA's values in y.")
  if(p >= 1 | p <= 0) stop("p must be a real number in (0,1).")
  if(abs(mu) ==Inf) stop("mu must be a finite real number.")
  
  dif = y - mu
  return((sum(p*dif[dif>0]) - sum((1-p)*dif[dif<0])))
}

mleALD = function(y,initial=NA)
{
  if(is.na(initial) == TRUE)
  {
    p = 0.5
    OPT = optimize(f = minfunc,interval = c(-10000,10000),y=y,p=p)
    mu = OPT$minimum
    sigma = (1/length(y))*OPT$objective
    
  }else
  {
    if(abs(initial[1]) ==Inf) stop("mu must be a finite real number.")
    mu = initial[1]
    if(initial[2] <= 0) stop("sigma must be a positive number.")
    sigma = initial[2]
    if(initial[3] >= 1 | initial[3] <= 0) stop("p must be a real number in (0,1).")
    p = initial[3]
  }

  teta0 = c(mu,sigma,p)
  criterio = 1
  iter = 0
  while(max(criterio) > 10^-5)
  {
    p = optimize(f = likALD,interval = c(0.001,0.999),y=y,mu=mu,sigma=sigma,loglik=T,maximum = T)$maximum
    OPT = optimize(f = minfunc,interval = c(-10000,10000),y=y,p=p)
    mu = OPT$minimum
    sigma = (1/length(y))*OPT$objective
    tetaT = c(mu,sigma,p)
    criterio = (teta0-tetaT)^2
    teta0 = tetaT
    iter = iter +1
  }
  
  return(list(iter = iter, par = teta0))
}

# y = rALD(10000,mu = -323,sigma = 40,p = 0.9)
# res = mleALD(y)
# seqq = seq(-3000,0)
# dens = dALD(y=seqq,mu=res$par[1],sigma=res$par[2],p=res$par[3])
# hist(y,breaks=50,freq = F,ylim=c(0,max(dens)))
# lines(seqq,dens,type="l",lwd=2,col="red",xlab="x",ylab="f(x)", main="ALD Density function")

Try the ald package in your browser

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

ald documentation built on April 5, 2021, 1:06 a.m.