R/dist.R

Defines functions distlist ed_skeleton ed_parse ft_stdt b_show fn_stdt fdist_stdt param_score_stdt stdnormmoment stdnormabsmoment stdtmoment stdtabsmoment tabsmoment fdist_stdnorm

Documented in b_show distlist ed_parse ed_skeleton fdist_stdnorm fdist_stdnorm fdist_stdt fn_stdt fn_stdt ft_stdt param_score_stdt stdnormabsmoment stdnormmoment stdtabsmoment stdtmoment tabsmoment

## Do not edit this file manually.
## It has been automatically generated from *.org sources.

fdist_stdnorm <- function(){
    param_flag <- FALSE
    f <- function(){
            list(  pdf = "dnorm"
                 , cdf = "pnorm"
                 , rand = "rnorm"
                 , logpdf = function(x) log(dnorm(x)) # this is lazy, of course
                 , Fscore = function(x) -x
                 , xFscore = function(x) -x^2
                 , Parscore = function() 0           #  todo: may be stop()?
                 , get_param = function() numeric()  #  todo: may be stop()?
                 , set_param = function() stop("Standard normal has no parameters to set.")
                 , any_param = function() param_flag
                 , show      = function()  "Standard normal distribution"
                 )
        }
    f()
}

dist_norm <- fdist_stdnorm()

# get_noise_dist <- function(x){
#     wrk <- noise_dists(x)
#     if(is.character(x)){
#         wrk <- noise_dists[x]   # should work for length(x)>1 too. (I think)
#
#     }
#
#
# }

## This chunk was in new.R
# from Mathematica

tabsmoment <- function(nu, k){         # exists if k < nu   E|x|^k za t-dist
    if(any(k >= nu)){
        res <- rep_len(Inf, length(k))
        ind <- k < nu 
        m <- k[ind]
        res[ind] <- nu^(m/2) * gamma((1 + m)/2) * gamma((nu-m)/2) /
                    (sqrt(pi) * gamma(nu/2))
        res
    }else
        nu^(k/2) * gamma((1 + k)/2) * gamma((nu-k)/2) / (sqrt(pi) * gamma(nu/2))
}

stdtabsmoment <- function(nu, k){   # nu>2 here      # exists if k < nu   E|x|^k za t-dist
    ## 2018-09-02 return Inf for k >= nu
    if(any(k >= nu)){
        res <- rep_len(Inf, length(k))
        ind <- k < nu 
        m <- k[ind]
        res[ind] <- (nu-2)^(m/2) * gamma((1 + m)/2) * gamma((nu-m)/2) /
                    (sqrt(pi) * gamma(nu/2))
        res
    }else
        (nu-2)^(k/2) * gamma((1 + k)/2) * gamma((nu-k)/2) / (sqrt(pi) * gamma(nu/2))
}
stdtmoment <- function(nu, k){   # nu>2 here      # exists if k < nu   E|x|^k za t-dist
    res <- rep_len(NaN, length(k))
    res[k == 0] <- 1                  # zero'th moment is = 1
    res[k %% 2 == 1] <- 0  # odd moments are zero (arbitrary, but principal value of the integral)
    i <- which(k %% 2 == 0  &  k > 0) # even moments are non-zero
    res[i] <- stdtabsmoment(nu, k[i])
    res
}

## for k > -1
stdnormabsmoment <- function(k)
    2^(k/2) * gamma((1 + k)/2) / sqrt(pi)

## this gives NaN's for fractional k's, TODO: make similar to stdt above?
stdnormmoment <- function(k)
    2^(k/2-1) * (1 + (-1)^k) * gamma((1 + k)/2) / sqrt(pi)

param_score_stdt <- function(x,nu){
   1/2*(1/(2-nu) + x^2*(1+nu)/((nu-2)*(x^2+nu-2)) - log(1+x^2/(nu-2))
           - digamma(nu/2) + digamma((nu+1)/2) )
}

fdist_stdt <- function(df, fixed = TRUE){# currently uses functions from fGarch; slagam
                                         # `fGarch' v "Depends"
    nu <- df    # todo: nu must be > 2 !
    param_flag <- !fixed
    f <- function(){
            list(  pdf = function(x) dstd(x, nu = nu)
                 , cdf = function(x) pstd(x, nu = nu)
                 , rand = function(n) rstd(n, nu = nu)
                 , logpdf = function(x) log(dstd(x, nu = nu))
                 , Fscore = function(x) - x*(1+nu)/(x^2+nu-2)
                 , xFscore = function(x) - x^2*(1+nu)/(x^2+nu-2)
                 , Parscore = function(x) param_score_stdt(x,nu)
                 , get_param = function() nu
                 , set_param = function(x) nu <<- x
                 , any_param = function() param_flag
                 , show = function() paste("Student t with", format(nu, digits = 4), "df")
                 )
        }
    f()
}

                                                              # todo: remove argument "fixed"?
                                  # note: df is a vector here!
fn_stdt <- function(df, fixed = TRUE){# currently uses functions from fGarch;
    nu <- df                                                          # todo: nu must be > 2 !
    param_flag <- !fixed                                           # todo: individual flags...
    res <- vector("list", length(df))
    creator_fun <- "fn_stdt"                   # todo: this may not be a long term solution...
    brief_show = function(){
        # cat("A list of the following distributions:\n")
        wrk <- sapply(res, function(x) x$show())
        cat(paste("\t", "Component ", seq_along(wrk), ": ", wrk, sep=""), sep="\n")
    }

    for(i in seq_along(nu)){
        local({
              k <- i
              res[[i]] <<- (function(){
                              list(  pdf = function(x) dstd(x, nu = nu[k])
                                   , cdf = function(x) pstd(x, nu = nu[k])
                                   , rand = function(n) rstd(n, nu = nu[k])
                                   , logpdf = function(x) log(dstd(x, nu = nu[k]))
                                   , Fscore = function(x) - x*(1+nu[k])/(x^2+nu[k]-2)
                                   , xFscore = function(x) - x^2*(1+nu[k])/(x^2+nu[k]-2)
                                   , Parscore = function(x) param_score_stdt(x,nu[k])
                                   , get_param = function() nu[k]
                                   , set_param = function(x) nu[k] <<- x
                                   , any_param = function() param_flag
                                   , show = function()
                                           paste("Student t with",
                                                 format(nu[[k]], digits = 4), "df")
                                   )
                          })()   # evaluated!; notice the parentheses around 'function'
                                 #   otherwise it is not evaluated (and no error is signaled)
          })}
    res
}

b_show <- function(x){
                                   # 2012-10-31 wrk <- x[[c(1,1)]] does not work if it is a
                                   #           string (e.g. "dnorm") or non-local function
    wrk <- x[[1]]$show  # todo: this may still fail; look for a better alternative

    if(is.function(wrk)){
        # f <- get("brief_show", envir = environment(wrk))
        f <- mget("brief_show", envir = environment(wrk), mode="function",
                  inherits = TRUE, ifnotfound = list(NULL))[[1]] # mget returns a list
        if(!is.null(f))
            f()
        else {
            lapply(seq_along(x), function(i) cat("\tComponent ", i, ": ", x[[i]]$show(), "\n", sep = ""))
            invisible(x)
        }
    ## }else if(is.list(x)){
    ##     lapply(x, function(d) d$show())
    ##     invisible(x)
    }else{
        show(x)
    }
}

# currently uses functions from fGarch;
                                  # note: df is a vector here!
ft_stdt <- function(df, fixed = TRUE, n = length(df), tr = function(x, k) x[k] ){
#cat("a: ", n, "\n")    
    ## 2021-06-26
    ##     was: if(is.null(n))  n <- if(missing(tr) || is.function(tr)) length(df) else length(tr)
    ##     but that leaves n = 0 when tr is NULL
    if(is.null(n))  n <- if(missing(tr) || is.function(tr) || is.null(tr)) length(df) else length(tr)
#cat("b: ", n, "\n")

    if(is.numeric(tr)){        # elem of 'tr' are indices in 'nu'
        nuind <- tr
        tr <- function(x,k) x[nuind[k]]
    }

    if(is.null(tr)) tr <- function(x, k) x[k]

    nu <- df                                                          # todo: nu must be > 2 !
    param_flag <- any(!fixed)                                      # todo: individual flags...
    param_flags <- !fixed            # 2012-10-25 individual flags; fixed may be vector here
    if(length(param_flags)==1)
        param_flags <- rep(param_flags, length(nu))

    get_non_fixed <- function() nu[param_flags]
    set_non_fixed <- function(x) nu[param_flags] <<- x

    get_nu        <- function() nu
    set_nu        <- function(x) nu <<- x

    creator_fun <- "ft_stdt"                   # todo: this may not be a long term solution...
    brief_show = function(){
        # cat("A list of the following distributions:\n")
        wrk <- sapply(res, function(x) x$show())
        cat(paste("\t", "Component ", seq_along(wrk), ": ", wrk, sep=""), sep="\n")
        # todo: put here info about common param and/or tr
    }

    res <- vector("list", n)
#cat("c: ", n, "\n")    
    for(i in 1:n){  # assumes n >= 1
        local({
              k <- i
              res[[i]] <<- (function(){
                              list(  pdf = function(x) dstd(x, nu = tr(nu,k))
                                   , cdf = function(x) pstd(x, nu = tr(nu,k))
                                   , rand = function(n) rstd(n, nu = tr(nu,k))
                                   , logpdf = function(x) log(dstd(x, nu = tr(nu,k)))
                                   , Fscore = function(x) - x*(1+tr(nu,k))/(x^2+tr(nu,k)-2)
                                   , xFscore = function(x) - x^2*(1+tr(nu,k))/(x^2+tr(nu,k)-2)
                                   , Parscore = function(x) param_score_stdt(x,tr(nu,k))
                                   , get_param = function() tr(nu,k)

                                    # not clear what to do here.
                                    # , set_param = function(x) nu[k] <<- x

                                   # na tova mozhe da se pridade smisal, no zasega macham.
                                   #      (taka nyama a ima opasnost ot confusion)
                                   ##
                                   ## 2021-06-26: reinstating but returning the overall flag
                                   ##             indicating if any dist in the distlist has
                                   ##             a parameter to estimate
                                   #        any_param = function() param_flags[k]
                                   , any_param = function() param_flag

                                   , show = function()
                                                   paste("Student t with",
                                                         format(tr(nu,k), digits = 4), "df")
                                   )
                          })()   # evaluated!; notice the parentheses around 'function'
                                 #   otherwise it is not evaluated (and no error is signaled)
          })}
    res
}

ed_stdnorm <-
    c(  pdf       = '"dnorm"'
      , cdf       = '"pnorm"'
      , rand      = '"rnorm"'
      , logpdf    = 'function(x) log(dnorm(x))'    # this is lazy, of course
      , Fscore    = 'function(x) -x'
      , xFscore   = 'function(x) -x^2'
      , Parscore  = 'function() 0'                         #  todo: may be stop()?
      , get_param = 'function() numeric()'     #  todo: may be stop()?
      , set_param = 'function() stop("Standard normal has no parameters to set.")'
      , any_param = 'function() param_flag'
      , show      = 'function()  "Standard normal distribution"'
      , moment    = 'function(x) stdnormmoment(x)'   # 2012-11-10 new
      )

ed_stdt0 <-
    c(  pdf       = 'function(x) dstd(x, nu = nu)'
      , cdf       = 'function(x) pstd(x, nu = nu)'
      , rand      = 'function(n) rstd(n, nu = nu)'
      , logpdf    = 'function(x) log(dstd(x, nu = nu))'
      , Fscore    = 'function(x) - x*(1+nu)/(x^2+nu-2)'
      , xFscore   = 'function(x) - x^2*(1+nu)/(x^2+nu-2)'
      , Parscore  = 'function(x) param_score_stdt(x,nu)'
      , get_param = 'function() nu'
      , set_param = 'function(x) nu <<- x'
      , any_param = 'function() param_flag'
      )

ed_stdt1 <-
    c(  pdf       = 'function(x) dstd(x, nu = nu[k])'
      , cdf       = 'function(x) pstd(x, nu = nu[k])'
      , rand      = 'function(n) rstd(n, nu = nu[k])'
      , logpdf    = 'function(x) log(dstd(x, nu = nu[k]))'
      , Fscore    = 'function(x) - x*(1+nu[k])/(x^2+nu[k]-2)'
      , xFscore   = 'function(x) - x^2*(1+nu[k])/(x^2+nu[k]-2)'
      , Parscore  = 'function(x) param_score_stdt(x,nu[k])'
      , get_param = 'function()  nu[k]'
      , set_param = 'function(x) nu[k] <<- x'
      , any_param = 'function()  param_flag'
      , show      = 'function()  paste("Student t with", format(nu[[k]], digits = 4), "df")'
      )

ed_stdt <-
    c(  pdf       = 'function(x) dstd(x, nu = tr(nu,k))'
      , cdf       = 'function(x) pstd(x, nu = tr(nu,k))'
      , rand      = 'function(n) rstd(n, nu = tr(nu,k))'
      , logpdf    = 'function(x) log(dstd(x, nu = tr(nu,k)))'
      , Fscore    = 'function(x) - x*(1+tr(nu,k))/(x^2+tr(nu,k)-2)'
      , xFscore   = 'function(x) - x^2*(1+tr(nu,k))/(x^2+tr(nu,k)-2)'
      , Parscore  = 'function(x) param_score_stdt(x,tr(nu,k))'
      , get_param = 'function() tr(nu,k)'
      , show      = 'function() paste("Student t with", format(tr(nu,k), digits = 4), "df")'
      , moment    = 'function(x) stdtmoment(tr(nu,k), x)'   # 2012-11-10 new
      ## 2021-08-23: add the below any_param(). this is misleading but it is the same in the 
      ##             N(0,1), which doesn't have param. It seems that the intent had been this
      ##              to be a flag about any param in any component.
      ##      TODO: this needs big consolidation. In particular there should be distinction 
      ##            between any_param for individual component and any_param for all 
      ##            (in that case probably 'parflag' could be an attribute.
      ##     But if any_param() is used consistently in the current sense, this may be ok.
      , any_param = 'function()  param_flag'
      )
## These comments are from 2012 or older:
      # , 'not clear what to do here.'
      # , 'set_param = function(x) nu[k] <<- x'
      # ', '  ## 2017-04-21 this was #', ', causing devtools to create documentation.
      # , 'na tova mozhe da se pridade smisal, no zasega macham.'
      # , '     (taka nyama a ima opasnost ot confusion)'
      # , 'any_param = function() param_flags[k]'

ed_parse <- function(s){
    nams <- allNames(s)
    fillers <- ifelse(nams == "", "", "=")
    wrk <- paste(nams, fillers, s, collapse = ",\n")

    text <- c("function(){", "list(", wrk, ")", "}")
    res <- parse(text = text)
    res
}

ed_src <- list(  stdnorm = ed_parse(ed_stdnorm)
               , stdt    = ed_parse(ed_stdt)
               )

ed_nparam <- c(stdnorm = 0, stdt = 1)

ed_skeleton <- function(df, fixed = FALSE, n = length(df), tr = NULL){
    if(is.null(tr)) tr <- function(x, k) x[k]
    else if(is.numeric(tr)){        # elem of 'tr' are indices in 'nu'
        nuind <- tr
        tr <- function(x,k) x[nuind[k]]
    }

    if(is.null(n))
        n <- if(is.function(tr)) length(df) else length(tr)

    nu <- df
    if(length(nu) == 0){  # arg. 'fixed' is irrelevant in this case
        param_flag  <- logical(0)
        param_flags <- logical(0)
    }else{
        param_flag <- any(!fixed)
        param_flags <- !fixed          # 2012-10-25 individual flags; fixed may be vector here
    }
    if(length(param_flags)==1)
        param_flags <- rep(param_flags, length(nu))

    get_non_fixed <- function() nu[param_flags]
    set_non_fixed <- function(x) nu[param_flags] <<- x

    get_nu        <- function() nu
    set_nu        <- function(x) nu <<- x

    creator_fun <- "ed_skeleton"               # todo: this may not be a long term solution...
    brief_show = function(){
        # cat("A list of the following distributions:\n")
        wrk <- sapply(res, function(x) x$show())
        cat(paste("\t", "Component ", seq_along(wrk), ": ", wrk, sep=""), sep="\n")
        # todo: put here info about common param and/or tr
    }

    res <- vector("list", n)
    for(i in 1:n){  # assumes n >= 1
        local({
              k <- i
              res[[k]] <<- function(x){
                               res[[k]] <<- (eval(x,parent.env(environment())))()
                               res[[k]] }
          })}
    res
}

                                        # 2012-10-28 new
distlist <- function(type, param, ncomp = NULL, fixed = FALSE, tr = NULL,  ...){
    if(is.function(type))
        return( list(generator = type, param = param) )

    trflag <- is.null(tr)

    if(is.character(type)){
        dist <-
        if(length(type) == 1){
            switch(type
                   , "stdnorm" =
                   , "stdnormal" =
                   , "stdGaussian" = list(generator = fdist_stdnorm, param = numeric(0))

                   , "stdt" = {
                                                 # interesting questions about evaluation here
                                      # copy the arguments locally to avoid (reduce) surprises
                       param <- if(missing(param) || is.null(param)) NA_integer_ else param
                       fixed <- fixed
                       if(length(param) == 1  && trflag)
                           function(par) fdist_stdt(par, fixed=fixed)
                       else{
                           ncomp <- ncomp
                           tr <- tr
                           function(par)
                               ft_stdt(par, fixed = fixed, n = ncomp, tr = tr)
                       }
                     }
                   ## default
                   , stop("unknown noise distribution.")
                   )

        }else{
                   ## tova pravi dist, no az iskam generator function.
                   ## dist <- ed_skeleton(param, fixed = fixed, n = length(type), tr = tr)
                   ## for(i in seq_along(type)){
                   ##     dist[[i]] <- dist[[i]](ed_src[[ type[i] ]])
                   ## }
            if(missing(param)){
                # param <- numeric(0)      #       params NA's would be better
                np <- sapply(type, function(x) ed_nparam[x])
                param <- rep(NA_real_, sum(np))
                if(is.null(tr)){
                    tr <- cumsum(np)# todo: this works only if the dist. have 0 or 1 param.
                }                   # but to allow more I need to check other places, as well.
            }                       # dist. with zero params also get an index here but this
                                    # is never referenced, so should not be a problem.
                                    # (note also that np[1]==0 is a special case of this case)

            function(par = numeric(0)){
                dist <- ed_skeleton(par, fixed = fixed, n = length(type), tr = tr)
                for(i in seq_along(type)){
                    dist[[i]] <- dist[[i]](ed_src[[ type[i] ]])
                }
                dist
            }
        }
    }else{
        stop("This branch of distlist() is not ready")
    }

    list(generator = dist, param = param)
}
GeoBosh/mixAR documentation built on May 9, 2022, 7:36 a.m.