FunctionMatrixCovariance
// 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 In probability theory and statistics, a covariance matrix (also known as auto-covariance matrix, dispersion matrix, variance matrix, or variance–covariance matrix) is a square matrix giving the covariance between each pair of elements of a given random vector.
// Intuitively, the covariance matrix generalizes the notion of variance to multiple dimensions. As an example, the variation in a collection of random points in two-dimensional space cannot be characterized fully by a single number, nor would the variances in the `x` and `y` directions contain all of the necessary information; a `2 × 2` matrix would be necessary to fully characterize the two-dimensional variation.
// Any covariance matrix is symmetric and positive semi-definite and its main diagonal contains variances (i.e., the covariance of each element with itself).
// The covariance matrix of a random vector `X` is typically denoted by `Kxx`, `Σ` or `S`.
// ~wikipedia.
library("FunctionMatrixCovariance")

// @function Estimate Covariance matrix with provided data. 
// @param M		`matrix<float>`		Matrix with vectors in column order.
// @returns Covariance matrix of provided vectors.
// 
// ---
// https://en.wikipedia.org/wiki/Covariance_matrix
// https://numpy.org/doc/stable/reference/generated/numpy.cov.html
export method cov (matrix<float> M, bool bias = false) =>
	int _rows = M.rows()
	int _cols = M.columns()
	switch 
		_rows <= 0 => runtime.error('Data vectors are empty.')
		_cols <  1 => runtime.error('Need atleast one data vector.')
	matrix<float> _S = matrix.new<float>(_cols, _cols, na)
	for _col = 0 to _cols - 1
		_Mcol = M.col(_col)
		for _i = _col to _cols - 1
			_c = _Mcol.covariance(M.col(_i), bias)
			_S.set(_col, _i, _c)
			_S.set(_i, _col, _c)
		_S.set(_col, _col, _Mcol.variance(bias))
	_S

// TEST 20230828 RS
// if barstate.islastconfirmedhistory
// 	iv = matrix.new<float>(3, 0, 0.0)
// 	iv.add_col(0, array.from(1.0, 0.5, 0.5))
// 	iv.add_col(1, array.from(0.5, 1.0, 0.5))
// 	iv.add_col(2, array.from(0.5, 0.5, 1.0))
// 	// [[1, 0.5, 0.5], [0.5, 1, 0.5], [0.5, 0.5, 1]]
// 	a = matrix.new<float>(5, 0, 0.0)
// 	a.add_col(0, array.from(1.0,2,3,4,5))
// 	a.add_col(1, array.from(6.0,7,8,9,40))
// 	b = matrix.new<float>(3, 0, 0.0)
// 	b.add_col(0, array.from(1.0,1,3))
// 	b.add_col(1, array.from(2.0,2,3))
// 	log.info('\n{0}', cov(iv))
// 	log.info('\n{0}', cov(a))
// 	log.info('\n{0}', cov(b))

