Probability
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © RicardoSantos

//@version=5

// @description: a library of functions to deal with probabilities.
library(title='Probability')

// imports:
import RicardoSantos/MathExtension/1 as me
// me.log2()
import RicardoSantos/ArrayGenerate/1 as agen
// agen.sequence_float()


//#region   Helper functions:

// @function Calculate the power of a square matrix
// @param M         Square matrix.
// @param power     Power to raise matrix by.
// @returns Raised matrix.
// 
// ---
// https://www.geeksforgeeks.org/finding-the-probability-of-a-state-at-a-given-time-in-a-markov-chain-set-2/
matrix_power (matrix<float> M, int power) =>
    int _dim = M.rows()
    if _dim != M.columns()
        runtime.error('Matrix is not square!')
    matrix<float> _A = matrix.new<float>(_dim, _dim, 0.0)
    matrix<float> _C = M.copy()
    // for _i = 1 to power
    //     _C := f_matrix_multiply_2d(_C, matrix_a, dimension, dimension, dimension, dimension)
    for _i = 0 to _dim - 1
        _A.set(_i, _i, 1.0)
    int _p = power
    for _i = 0 to 999999
        if _p > 0
            if (_p % 2) == 1
                _A := matrix.mult(_A, _C)
                // _A := f_matrix_multiply_2d(matrix_a=_A, matrix_b=_C, dimension_ax=dimension, dimension_ay=dimension, dimension_bx=dimension, dimension_by=dimension)
            _C := matrix.mult(_C, _C)
            // _C := f_matrix_multiply_2d(matrix_a=_C, matrix_b=_C, dimension_ax=dimension, dimension_ay=dimension, dimension_bx=dimension, dimension_by=dimension)
            _p /= 2
            int(na)
        else
            break
            int(na)
    _A


// @function replaces a labels text with the representation of a 2d matrix.
// @param label_id      Label to update.
// @param M             Data source.
// @returns Updates provided label text.
method label_matrix (label label_id, matrix<float> M)=>
    int _Mrows = M.rows()
    int _Mcols = M.columns()
    label.set_text(id=label_id, text=str.tostring(M, '#.###'))




//#endregion
//#region   Probability:
// reference:
//      https://gist.github.com/Daniel-Hug/f0fee168155f1af72c390030cc9be463
//      https://support.minitab.com/en-us/minitab-express/1/help-and-how-to/basic-statistics/probability-distributions/how-to/cumulative-distribution-function-cdf/methods-and-formulas/methods-and-formulas/
//      https://support.minitab.com/en-us/minitab-express/1/help-and-how-to/basic-statistics/probability-distributions/supporting-topics/basics/calculating-probabilities-for-different-distributions/

// @function Complementary error function
// @param value     Value to test.
// @returns float
// 
// ---
// `.erf(0.95) = 0.820891`
// 
// ---
// https://www.johndcook.com/blog/csharp_erf/
// https://reference.wolfram.com/language/ref/Erf.html
export erf(float value) =>
    // constants:
    float _a1 = 0.254829592,  float _a2 = -0.284496736, float _a3 = 1.421413741
    float _a4 = -1.453152027, float _a5 = 1.061405429,  float _p  = 0.3275911
    float _x = math.abs(value)
    float _t = 1.0 / (1.0 + _p * _x)
    // erf(-_x) = -erf(x)
    float _y = 1.0 - ((((_a5 * _t + _a4) * _t + _a3) * _t + _a2) * _t + _a1) * _t * math.exp(-_x * _x)
    math.sign(value) * _y


// @function Computes the inverse error function using the Mc Giles method, sacrifices accuracy for speed.
// @param value     value to test, expected value in (-1.0, 1.0) range.
// @returns float
// 
// ---
// `.ierf_mcgiles(0.820891) = 0.72096 // correct = 0.95` \
// `.ierf_mcgiles(0.6) = -0.17696 // correct = 0.595116`
// 
// ---
// http://www.mimirgames.com/articles/programming/approximations-of-the-inverse-error-function/
// https://keisan.casio.com/exec/system/1180573448
export ierf_mcgiles(float value) =>
    float _p = 0.
    float _w = -math.log((1.0 - value) * (1.0 + value))
    if _w < 0.5
        _w -= 2.5
        _p += ((((((((2.81022636e-08 * _w) + 3.43273939e-07 * _w) + -3.5233877e-06 * _w) + -4.39150654e-06 * _w) + 0.00021858087 * _w) + -0.00125372503 * _w) + -0.00417768164 * _w) + 0.246640727 + _w) + 1.50140941
    else
        _w := math.sqrt(nz(_w)) - 3.0
        _p += ((((((((-0.000200214257 * _w) + 0.000100950558 * _w) + 0.00134934322 * _w) + -0.00367342844 * _w) + 0.00573950773 * _w) + -0.0076224613 * _w) + 0.00943887047 * _w) + 1.00167406 * _w) + 2.83297682
    _p * value


// @function computes the inverse error function using the Newton method with double refinement.
// @param value     value to test, expected value in (-1.0, 1.0) range.
// @returns float
// 
// ---
//  http://www.mimirgames.com/articles/programming/approximations-of-the-inverse-error-function/
// 
// ---
// `.ierf(0.820891) = na` //https://reference.wolfram.com/language/ref/InverseErf.html \
// `.ierf_mcgiles(0.6) = na` // correct = 0.595116 // https://keisan.casio.com/exec/system/1180573448
export ierf_double(float value) =>
    float _r = 0.
    float _z = math.sign(value) * value
 
    if _z <= 0.7
        float _x2 = _z * _z
        _r += _z * (((-0.140543331 * _x2 + 0.914624893) * _x2 - 1.645349621) * _x2 + 0.886226899)
        _r /= (((0.012229801 * _x2 - 0.329097515) * _x2 + 1.442710462) * _x2 - 2.118377725) * _x2 + 1.
    else
        float _y = math.sqrt( -math.log((1 - _z) / 2))
        _r += (((1.641345311 * _y + 3.429567803) * _y - 1.62490649) * _y - 1.970840454)
        _r /= ((1.637067800 * _y + 3.543889200) * _y + 1.)
    _r *= math.sign(value)
    _z *= math.sign(value)
    
    _r -= (erf(value=_r) - _z) / (2 / math.sqrt(math.pi) * math.exp(-_r * _r))
    _r -= (erf(value=_r) - _z) / (2 / math.sqrt(math.pi) * math.exp(-_r * _r)) //Comment out if you do not want single refinement


// @function computes the inverse error function using the Newton method.
// @param value     value to test, expected value in (-1.0, 1.0) range.
// @returns float
// 
// ---
//  http://www.mimirgames.com/articles/programming/approximations-of-the-inverse-error-function/
// 
// ---
//      .ierf(0.820891) =  0.95 //https://reference.wolfram.com/language/ref/InverseErf.html
//      .ierf_mcgiles(0.6) = 0.59512 // correct = 0.595116 // https://keisan.casio.com/exec/system/1180573448
export ierf(float value) =>
    float _r = 0.
    float _z = math.sign(value) * value
 
    if _z <= 0.7
        float _x2 = _z * _z
        _r += _z * (((-0.140543331 * _x2 + 0.914624893) * _x2 - 1.645349621) * _x2 + 0.886226899)
        _r /= (((0.012229801 * _x2 - 0.329097515) * _x2 + 1.442710462) * _x2 - 2.118377725) * _x2 + 1.
    else
        float _y = math.sqrt(-math.log((1 - _z) / 2))
        _r += (((1.641345311 * _y + 3.429567803) * _y - 1.62490649) * _y - 1.970840454)
        _r /= ((1.637067800 * _y + 3.543889200) * _y + 1.)
    _r *= math.sign(value)
    _z *= math.sign(value)
    
    _r -= (erf(_r) - _z) / (2 / math.sqrt(math.pi) * math.exp(-_r * _r))
    // _r -= (erf(_r) - _z) / (2 / math.sqrt(math.pi) * math.exp(-_r * _r)) //Comment out if you do not want single refinement


// @function probability that the event will not occur.
// @param probability       Probability of event, expected value in (0.0, 1.0) range.
// @returns float
// 
// ---
// `probability_complement(0.25) = 0.75` \
// `probability_complement(0.67) = 0.23`
export complement(float probability) =>
    1.0 - probability


// @function Gini Inbalance or Gini index for a given probability.
// @param probability       Probability of event, expected value in (0.0, 1.0) range.
// @returns float
// 
// ---
// https://victorzhou.com/blog/gini-impurity/
// https://www.tradingview.com/script/IlAZiO7E-Function-Entropy-Gini-Index/
export entropy_gini_impurity_single(float probability) =>
    probability * complement(probability=probability)


// @function Gini Inbalance or Gini index for a series of events.
// @param events        List with event probability's, expected values to be in (0.0, 1.0) range.
// @returns float
// 
// ---
// https://victorzhou.com/blog/gini-impurity/
// https://www.tradingview.com/script/IlAZiO7E-Function-Entropy-Gini-Index/
export entropy_gini_impurity(float[] events) =>
    int _size = array.size(id=events)
    float _sum = 0.
    if _size > 0
        for _i = 0 to _size - 1
            _p = array.get(id=events, index=_i)
    	    _sum += entropy_gini_impurity_single(probability=_p)
	_sum


// @function Entropy information value of the probability of a single event.
// @param probability       Probability of event, expected value in (0.0, 1.0) range.
// @returns float, value as bits of information.
// 
// ---
// https://machinelearningmastery.com/what-is-information-entropy/
// https://www.tradingview.com/script/JFXtLbSo-Function-Shannon-Entropy/
export entropy_shannon_single(float probability) =>
    -me.log2(value=probability + 1.0e-15)


// @function Entropy information value of a distribution of events.
// @param events        List with event probability's, expected values to be in (0.0, 1.0) range.
// @returns float
// 
// ---
// https://machinelearningmastery.com/what-is-information-entropy/
// https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.entropy.html
// https://www.tradingview.com/script/JFXtLbSo-Function-Shannon-Entropy/
export entropy_shannon(float[] events) =>
    int _size = array.size(id=events)
    float _sum = 0.
    if _size > 0
        for _i = 0 to _size - 1
            _p = array.get(id=events, index=_i)
    	    _sum += _p * entropy_shannon_single(probability=_p)
	_sum


// @function Calculates Chebyshev Inequality.
// @param n_stdeviations    Value expected to be positive, over or equal to 1.0.
// @returns float
// 
// ---
// https://www.tradingview.com/script/ndG6Dffb-Function-Probability-Chebyshev-Inequality/
export inequality_chebyshev(float n_stdeviations=1.0) =>
    n_stdeviations < 1 ? 0 : 1. - 1. / math.pow(base=n_stdeviations, exponent=2.0)


// @function Calculates Chebyshev Inequality.
// @param mean      Mean of a distribution.
// @param std       standard deviation of a distribution.
// @returns float
// 
// ---
// https://www.tradingview.com/script/ndG6Dffb-Function-Probability-Chebyshev-Inequality/
export inequality_chebyshev_distribution(float mean, float std) =>
    inequality_chebyshev(n_stdeviations=mean / std)


// @function Calculates Chebyshev Inequality for a array of values.
// @param data_sample       List of numeric values.
// @returns float
// 
// ---
// https://www.tradingview.com/script/ndG6Dffb-Function-Probability-Chebyshev-Inequality/
export inequality_chebyshev_sample(float[] data_sample) =>
    float _mean = array.avg(id=data_sample)
    float _std = array.stdev(id=data_sample)
    inequality_chebyshev(n_stdeviations=_mean / _std)


// @function Probability that all arguments will happen when neither outcome
//  is affected by the other (accepts 1 or more arguments)
// @param events    List of event probabilities.
// @returns Probability.
// 
// ---
// `.intersection_of_independent_events([0.25, 0.25]) = 0.0625` \
// `.intersection_of_independent_events([0.3, 0.3, 0.3]) = 0.027`
export intersection_of_independent_events(float[] events) =>
    int _size = array.size(id=events)
    if _size > 0
        float _ret = array.get(id=events, index=0)
        if _size > 1
        	for _i = 1 to _size - 1
        		_ret *= array.get(id=events, index=_i)
    	_ret


// @function Probability that either one of the arguments will happen when neither outcome
//  is affected by the other (accepts 1 or more arguments)
// @param events    List of event probabilities.
// @returns Probability.
// 
// ---
// `.union_of_independent_events([0.25, 0.25]) = 0.4375`   \
// `.union_of_independent_events([0.3, 0.3, 0.3]) = 0.657` 
export union_of_independent_events(float[] events) =>
    int _size = array.size(id=events)
    if _size > 0
        float _ret = array.get(id=events, index=0)
        if _size > 1
        	for _i = 1 to _size - 1
        	    float _Ai = array.get(id=events, index=_i)
        		_ret += _Ai - nz(_ret) * _Ai
        		//arguments[i] - ret * arguments[i];
    	_ret


// @function Probabilities for each bin in the range of sample.
// @param sample    Samples to pool probabilities.
// @param n_bins    Number of bins to split the range
// @return Probability distribution.
// 
// ---
// `.mass_function([1, 2, 3, 4], 4) = [0.25, 0.25, 0.25, 0.25]`
export mass_function(float[] sample, int n_bins) =>
    float[] _probabilities = array.new_float(size=n_bins, initial_value=0.0)
    int _size = array.size(id=sample), float _fraction = 1.0 / _size
    float _min = array.min(id=sample), float _max = array.max(id=sample), float _range = _max - _min
    if _size > 0
        for _i = 0 to _size - 1
            int _expected_index = math.floor(((array.get(id=sample, index=_i) - _min) / _range) * (n_bins - 1))
            array.set(id=_probabilities, index=_expected_index, value=array.get(id=_probabilities, index=_expected_index) + _fraction)
    _probabilities


// @function Use the CDF to determine the probability that a random observation 
//  that is taken from the population will be less than or equal to a certain value.
//  Or returns the area of probability for a known value in a normal distribution.
// @param mean      Mean of the distribution.
// @param stdev     Standard deviation of the distribution.
// @param value     Limit value at which to stop.
// @returns Probability.
// 
// ---
// `.cumulative_distribution_function(5, 30, 25) = 0.747507...` \
// `.cumulative_distribution_function(30, 25, 1.4241) = 0.126512`
// 
// ---
// https://en.wikipedia.org/wiki/Cumulative_distribution_function
// https://stackoverflow.com/questions/5259421/cumulative-distribution-function-in-javascript
// https://www.wolframalpha.com/input/?i=CDF%28normal+distribution+with+mean%3D+5+and+sd+%3D+30%2C+25%29
export cumulative_distribution_function(float mean, float stdev, float value) =>
    _z = (value - mean) / math.sqrt(2.0 * math.pow(stdev, 2.0))
    0.5 * (1.0 + erf(_z))


// @function Transition matrix for the suplied distribution.
// @param distribution      Probability distribution. ex:. [0.25, 0.50, 0.25]
// @returns Transition matrix.
export transition_matrix(float[] distribution) =>
    int _size = array.size(id=distribution)
    if _size < 0
        runtime.error('Can not generate transition matrix from a empty distribution!')
    matrix<float> _T = matrix.new<float>(_size, _size, na)
    for _m = 0 to _size - 1
        for _n = 0 to _size - 1
            _T.set(_m, _n, distribution.get(_m) * distribution.get(_n))
    _T


// @function Probability of reaching target_state at target_step after starting from start_state
// @param transition_matrix     Probability transition matrix.
// @param dimension             Size of the matrix dimension.
// @param target_step           Number of steps to find probability.
// @returns Probability.
//
// ---
// https://www.geeksforgeeks.org/finding-the-probability-of-a-state-at-a-given-time-in-a-markov-chain-set-2/
export diffusion_matrix (matrix<float> transition_matrix, int target_step) =>
    matrix_power(transition_matrix, target_step)


// @function Probability of reaching target_state at target_step after starting from start_state
// @param transition_matrix     Probability transition matrix.
// @param start_state           State at which to start.
// @param target_state          State to find probability.
// @param target_step           Number of steps to find probability.
// @returns Probability.
//
// ---
// https://www.geeksforgeeks.org/finding-the-probability-of-a-state-at-a-given-time-in-a-markov-chain-set-2/
export state_at_time (matrix<float> transition_matrix, int start_state, int target_state, int target_step) =>
    matrix<float> _T = matrix_power(transition_matrix, target_step)
    _T.get(target_state - 1, start_state - 1)

//#endregion






// import RicardoSantos/Probability/1

var label test = label.new(bar_index, 0.0, '')
if barstate.islastconfirmedhistory
    t0 = erf(0.6)
    t1 = ierf(0.6)
    t2 = intersection_of_independent_events(array.new_float(3, 0.3))
    t3 = union_of_independent_events(array.new_float(3, 0.3))
    t4 = cumulative_distribution_function(30, 25, 1.4241)
    seq = agen.sequence_float(0.0, 1.0, 0.25) // [0.0, 0.25, 0.5, 0.75, 1.0]
    tm = transition_matrix(seq)
    label.set_xy(test, bar_index, 0.0)
    //label.set_text(test, tostring(t4))
    // label_matrix(test, transition_matrix(array.new_float(4, 0.25)))
    // label_matrix(test, transition_matrix(seq))
    //f_label_matrix_2d(test, f_matrix_multiply_2d(probability_transition_matrix(array.new_float(4, 0.25)), probability_transition_matrix(array.new_float(4, 0.25)), 4, 4, 4, 4), 4, 4)
    label_matrix(test, diffusion_matrix(tm, 25))
    // label.set_text(test, str.tostring(state_at_time(tm, 2, 3, 25)))

    // 20230830 RS: previous version had wrong dimension value (4 -> 5)..
    // log.info('\n{0}', Probability.state_at_time(Probability.transition_matrix(seq), 4, 2, 3, 25))
    // log.info('\n{0}\n\n{1}',
    //  diffusion_matrix(tm, 25), 
    //  Probability.diffusion_matrix(Probability.transition_matrix(seq), 5, 25))

