MathProbabilityDistribution
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © RicardoSantos
// Special Apreciation for help to: Anton Berlin, LucF.

//@version=5

// @description Probability Distribution Functions.
library(title='MathProbabilityDistribution')

// reference:
//      https://r-forge.r-project.org/scm/viewvc.php/pkg/fGarch/src/dist.f?view=markup&revision=5924&root=rmetrics
//      https://itl.nist.gov/div898/handbook/eda/section3/eda3661.htm
//      https://rdrr.io/cran/fGarch/src/R/dist-sstd.R
//      https://rdrr.io/rforge/fGarch/

import RicardoSantos/MathSpecialFunctionsGamma/1 as msf
// .gamma()
import RicardoSantos/Probability/1 as prob
// .erf() .cumulative_distribution_function()

// Common Constant Variables:
var float __2PI__   = 2.0 * math.pi
var float __SQ2PI__ = math.sqrt(__2PI__)

// Shared variables and utility:
int   __CYCLE4__  = bar_index % 81
float __CYCLE4F__ = (-40 + __CYCLE4__) * 0.1
plot(math.abs(__CYCLE4F__) == 0 ? 0 : na, title='Cycle 1 Marker', color=color.white, linewidth=2, style=plot.style_cross)
plot(math.abs(__CYCLE4F__) == 1 ? 0 : na, title='Cycle 1 Marker', color=color.maroon, linewidth=2, style=plot.style_cross)
plot(math.abs(__CYCLE4F__) == 2 ? 0 : na, title='Cycle 2 Marker', color=color.red, linewidth=2, style=plot.style_cross)
plot(math.abs(__CYCLE4F__) == 3 ? 0 : na, title='Cycle 3 Marker', color=color.orange, linewidth=3, style=plot.style_cross)
plot(math.abs(__CYCLE4F__) == 4 ? 0 : na, title='Cycle 4 Marker', color=color.yellow, linewidth=3, style=plot.style_cross)

// @function    Indexed names helper function.
// @param       idx    int, position in the range (0, 6).
// @returns     string, distribution name.
// usage:
//      .name(1)
// Notes:
// (0) => 'StdNormal'     
// (1) => 'Normal'        
// (2) => 'Skew Normal'   
// (3) => 'Student T'     
// (4) => 'Skew Student T'
// (5) => 'GED'           
// (6) => 'Skew GED'
export name (int idx) => //{
    switch idx
        (0) => 'StdNormal'     
        (1) => 'Normal'        
        (2) => 'Skew Normal'   
        (3) => 'Student T'     
        (4) => 'Skew Student T'
        (5) => 'GED'           
        (6) => 'Skew GED'
        =>     'StdNormal'
//}

// @function    Z-score helper function for x calculation.
// @param       position    float, position.
// @param       mean        float, mean.
// @param       deviation   float, standard deviation.
// @returns     float, z-score.
// usage:
//      .zscore(1.5, 2.0, 1.0)
export zscore (float position, float mean, float deviation) => //{
    (position - mean) / deviation
// test:
// plot(zscore(__CYCLE4F__, 0.0, 1.0)) // OK
//}

// @function    Standard Normal Distribution.
// @param       position    float, position.
// @returns     float, probability density.
// usage:
//      .std_normal(0.6)
export std_normal (float position) => //{
    math.exp(-math.pow(position, 2.0) / 2.0) / __SQ2PI__
// test:
// plot(std_normal(__CYCLE4F__)) // OK
//}

// @function    Normal Distribution.
// @param       position    float, position in the distribution.
// @param       mean        float, mean of the distribution, default=0.0 for standard distribution.
// @param       scale       float, scale of the distribution, default=1.0 for standard distribution.
// @returns     float, probability density.
// usage:
//      .normal(0.6)
export normal (float position, float mean=0.0, float scale=1.0) => //{
    math.exp(-math.pow(position - mean, 2.0) / (2.0 * math.pow(scale, 2.0))) / (scale * __SQ2PI__)
// test:
// plot(normal(__CYCLE4F__)) // OK
//}

// @function    Skew Normal Distribution.
// @param       position    float, position in the distribution.
// @param       skew        float, skewness of the distribution.
// @param       mean        float, mean of the distribution, default=0.0 for standard distribution.
// @param       scale       float, scale of the distribution, default=1.0 for standard distribution.
// @returns     float, probability density.
// usage:
//      .skew_normal(0.8, -2.0)
export skew_normal (float position, float skew, float mean=0.0, float scale=1.0) => //{
    2.0 * std_normal(position) * prob.cumulative_distribution_function(mean, scale, skew * position)
// test:
// plot(skew_normal(__CYCLE4F__, 0.5)) // OK
//}
    
// @function    Generalized Error Distribution.
// @param       position    float, position.
// @param       shape       float, shape.
// @param       mean        float, mean, default=0.0 for standard distribution.
// @param       scale       float, scale, default=1.0 for standard distribution.
// @returns     float, probability.
// usage:
//      .ged(0.8, -2.0)
export ged (float position, float shape, float mean=0.0, float scale=1.0) => //{
    float _1shape = 1.0 / shape
    float _g1shape = msf.Gamma(_1shape)
    float _g3shape = msf.Gamma(3.0 / shape)
    float _z = zscore(position, mean, scale)
    float _lambda = math.sqrt(math.pow(2.0, (-2.0 / shape)) * _g1shape / _g3shape)
    float _g = shape / (_lambda * (math.pow(2.0, 1.0 + _1shape) * _g1shape))
    _g * math.exp(-0.5 * math.pow(math.abs(_z / _lambda), shape)) / scale
// test:
//  run: https://rdrr.io/snippets/
//      library(fGarch)
//      x = seq(-4, 4, by=0.1)
//      #            x, m, s, n
//      plot(x, dged(x, 0, 1, 0.7))
//      box()
//      lines(x, dged(x, 0, 1, 1))
//      lines(x, dged(x, 0, 1, 2))
//      lines(x, dged(x, 0, 1, 100))
//      lines(x, dged(x, 0, 1, 2))
//      lines(x, dged(x, 0, 2, 2))
//      lines(x, dged(x, 0, 3, 2))
//      lines(x, dged(x, 0, 4, 2))
// plot(ged(__CYCLE4F__, 0.7))      // OK
// plot(ged(__CYCLE4F__, 1))        // OK
// plot(ged(__CYCLE4F__, 2))        // OK
// plot(ged(__CYCLE4F__, 100))      // OK
// plot(ged(__CYCLE4F__, 2, 0, 1))  // OK
// plot(ged(__CYCLE4F__, 2, 0, 2))  // OK
// plot(ged(__CYCLE4F__, 2, 0, 3))  // OK
// plot(ged(__CYCLE4F__, 2, 0, 4))  // OK
//}

// @function    Skew Generalized Error Distribution.
// @param       position    float, position.
// @param       shape       float, shape.
// @param       skew        float, skew.
// @param       mean        float, mean, default=0.0 for standard distribution.
// @param       scale       float, scale, default=1.0 for standard distribution.
// @returns     float, probability.
// usage:
//      .skew_ged(0.8, 2.0, 1.0)
export skew_ged (float position, float shape, float skew, float mean=0.0, float scale=1.0) => //{
    // common
    float _1shape  = 1.0 / shape
    float _2shape  = 2.0 / shape
    float _1skew   = 1.0 / skew
    float _g1shape = msf.Gamma(_1shape)
    float _g2shape = msf.Gamma(_2shape)
    float _g3shape = msf.Gamma(3.0 / shape)
    // standardize:
    float _lambda  = math.sqrt(math.pow(2.0, -_2shape) * _g1shape / _g3shape)
    // float _g = shape / (_lambda * math.pow(2.0, 1.0 + _1shape) * _g1shape) // not used
    float _m1 = math.pow(2.0, _1shape) * _lambda * _g2shape / _g1shape
    float _mu = _m1 * (skew - _1skew)
    float _xi2 = math.pow(skew, 2.0)
    float _m12 = math.pow(_m1, 2.0)
    float _sigma = math.sqrt((1.0 - _m12) * (_xi2 + 1.0 / _xi2) + 2.0 * _m12 - 1.0)
    float _z = zscore(position, mean, scale) * _sigma + _mu
    // compute:
    float _xxi = math.pow(skew, math.sign(_z))
    if (_z == 0.0)
        _xxi := math.pow(skew, 0.0)//math.sign(skew)//skew**0.0
    float _g = (2.0 / (skew + _1skew))
    float _density = _g * ged(_z / _xxi, shape)
    (_density * _sigma) / scale
// test:
//  run: https://rdrr.io/snippets/
//      library(fGarch)
//      x = seq(-4, 4, by=0.1)
//      #             x, m, s, n, xi
//      plot(x, dsged(x, 0, 1, 2, 1.5))
//      box()
//      lines(x, dsged(x, 0, 2, 2, 1.5))
//      lines(x, dsged(x, 0, 3, 2, 1.5))
//      lines(x, dsged(x, 0, 4, 2, 1.5))
// plot(skew_ged(__CYCLE4F__, 2, 1.5, 0, 1)) // OK
// plot(skew_ged(__CYCLE4F__, 2, 1.5, 0, 2)) // OK
// plot(skew_ged(__CYCLE4F__, 2, 1.5, 0, 3)) // OK
// plot(skew_ged(__CYCLE4F__, 2, 1.5, 0, 4)) // OK
//}

// @function Student-T Distribution.
// @param position float, position.
// @param shape float, shape.
// @param mean float, mean, default=0.0 for standard distribution.
// @param scale float, scale, default=1.0 for standard distribution.
// @returns float, probability.
// usage:
//      .student_t(0.8, 2.0, 1.0)
export student_t (float position, float shape, float mean=0.0, float scale=1.0) => //{
    float _z = zscore(position, mean, scale)
    float _s = math.sqrt(shape / (shape - 2))
    float _shape2 = (shape + 1.0) / 2.0
    float _gshape2 = msf.Gamma(_shape2)
    float _a = _gshape2 / math.sqrt(shape * math.pi)
    float _b = msf.Gamma(shape / 2.0) * math.pow((1.0 + math.pow(_z * _s, 2.0) / shape), _shape2)
    (_a / _b) * _s / scale
// test:
//  run: https://rdrr.io/snippets/
//      library(fGarch)
//      x = seq(-4, 4, by=0.1)
//      #            x, m, s, n
//      plot(x, dstd(x, 0, 1, 5))
//      box()
//      lines(x, dstd(x, 0, 2, 5))
//      lines(x, dstd(x, 0, 3, 5))
//      lines(x, dstd(x, 0, 4, 5))
// plot(student_t(__CYCLE4F__, 5, 0, 1)) // OK
// plot(student_t(__CYCLE4F__, 5, 0, 2)) // OK
// plot(student_t(__CYCLE4F__, 5, 0, 3)) // OK
// plot(student_t(__CYCLE4F__, 5, 0, 4)) // OK
//}

// @function    Skew Student-T  Distribution.
// @param       position    float, position.
// @param       shape       float, shape.
// @param       skew        float, skew.
// @param       mean        float, mean, default=0.0 for standard distribution.
// @param       scale       float, scale, default=1.0 for standard distribution.
// @returns     float, probability.
// usage:
//      .skew_student_t(0.8, 2.0, 1.0)
export skew_student_t (float position, float shape, float skew, float mean=0.0, float scale=1.0) => //{
    float _1skew = 1.0 / skew
    float _a = 0.5
    float _b = shape / 2.0
    float _beta = (msf.Gamma(_a) / msf.Gamma(_a + _b)) * msf.Gamma(_b)
    float _m1 = 2.0 * math.sqrt(shape - 2.0) / (shape - 1.0) / _beta
    float _mu = _m1 * (skew - _1skew)
    float _xi2 = math.pow(skew, 2.0)
    float _m12 = math.pow(_m1, 2.0)
    float _sigma = math.sqrt((1.0 - _m12) * (_xi2 + 1.0 / _xi2) + 2.0 * _m12 - 1.0)
    float _z = position * _sigma + _mu
    float _xxi = math.pow(skew, math.sign(_z))
    if (_z == 0.0)
        _xxi := math.pow(skew, 0.0)//math.sign(skew)//skew**0.0
    float _g = 2.0 / (skew + _1skew)
    _g * student_t(_z / _xxi, shape, mean, scale) * _sigma 
// test:
//  run: https://rdrr.io/snippets/
//      library(fGarch)
//      x = seq(-4, 4, by=0.1)
//      #             x, m,    s, n, xi
//      plot(x, dsstd(x, 0, 0.25, 5, 1))
//      box()
//      lines(x, dsstd(x, 0, 1, 5, 2))
//      lines(x, dsstd(x, 0, 1, 5, 3))
// plot(skew_student_t(__CYCLE4F__, 5, 1, 0, 0.25)) // OK
// plot(skew_student_t(__CYCLE4F__, 5, 2, 0, 1)) // OK
// plot(skew_student_t(__CYCLE4F__, 5, 3, 0, 1)) // OK
//}

// @function Conditional Distribution.
// @param       distribution    string, distribution name.
// @param       position        float, position.
// @param       mean            float, mean, default=0.0 for standard distribution.
// @param       scale           float, scale, default=1.0 for standard distribution.
// @param       shape           float, shape.
// @param       skew            float, skew.
// @param       log             bool, if true apply log() to the result.
// @returns     float, probability.
// usage:
//      .select('StdNormal', __CYCLE4F__, log=true)
export select (string distribution='StdNormal', float position, float mean=0.0, float scale=1.0, float skew=na, float shape=na, bool log=false) => //{
    float _p = switch (distribution)
        ('StdNormal'        ) => std_normal(position / scale) / scale
        ('Normal'           ) => normal(position, mean, scale)
        ('Skew Normal'      ) => skew_normal(position, skew, mean, scale)
        ('Student T'        ) => student_t(position, shape, mean, scale)
        ('Skew Student T'   ) => skew_student_t(position, shape, skew, mean, scale)
        ('GED'              ) => ged(position, shape, mean, scale)
        ('Skew GED'         ) => skew_ged(position, shape, skew, mean, scale)
        => float(na)
    log ? math.log(_p) : _p
// test:
plot(select('StdNormal', __CYCLE4F__, log=true)) // OK
//}

