Cet article présente l'algorithme de Donald Knuth, qui est un bon algorithme de calcul de dispersion. Je me fie généralement à la bibliothèque pour le calcul de la distribution, et je n'ai pas considéré l'algorithme en particulier, donc j'espère pouvoir profiter de cette occasion pour approfondir mes connaissances.
Tout d'abord, je vais vous expliquer pourquoi vous avez besoin d'un tel algorithme. La distribution non biaisée est exprimée par la formule suivante.
\sigma^2 = \frac{1}{n(n-1)} \biggl(\sum_{i=1}^{n}x_i^2-\Bigl(\sum_{i=1}^{n}x_i\Bigr)^2\biggr) =\frac{1}{n-1}\Bigl(E[X^2]-E[X]^2\Bigr)
Maintenant, voici deux sommes cumulatives. 1. Somme cumulée de $ x $ lui-même et de 2. $ x ^ 2 $. Il y a actuellement deux problèmes principaux. Premièrement, la précision des calculs numériques n'est pas garantie lorsque $ x $ est un très grand nombre. Deuxièmement, nous devons refaire de gros calculs lorsque nous avons plus d'échantillons. Pour résoudre ces défis, examinons l'algorithme de Donald Knuth dans la section suivante.
Pour résoudre les problèmes mentionnés ci-dessus, cet algorithme a été conçu pour calculer la variance séquentiellement. L'algorithme est illustré ci-dessous. $ M_k $ représente la moyenne et $ s_k $ est le numérateur de la dispersion.
Initialisation: $ M_1 = x_1, S_1 = 0 $ Calculez l'équation graduelle suivante en $ k \ geq 2 $
\begin{align}
M_k &= M_{k-1} + \frac{(x_k - M_{k-1})}{k}\\
S_k &= S_{k-1} + (x_k - M_{k-1})*(x_k - M_k)
\end{align}
La variance non biaisée estimée à calculer est $ s ^ 2 = S_k / (k-1) $ à $ k \ geq 2 $.
Implémentons-le en Python pour nous assurer que l'algorithme précédent fonctionne réellement. Puisqu'il s'agit d'une expression graduelle, nous allons l'implémenter avec une fonction récursive. J'ai rendu possible la vérification de chaque instruction d'impression afin que vous puissiez voir comment elle est calculée séquentiellement.
def calc_var(x, k=0, M=0, s=0):
print('k=', k)
print('M=', M)
print('s=', s)
print('-----------------------')
if k == 0:
M = x[0]
s = 0
delta = x[k] - M
M += delta/(k+1)
s += delta*(x[k]-M)
k += 1
if k == len(x):
return M, s/(k-1)
else:
return calc_var(x, k=k, M=M, s=s)
x = [3.3, 5, 7.2, 12, 4, 6, 10.3]
print(calc_var(x))
Résultat d'exécution
k= 0
M= 0
s= 0
-----------------------
k= 1
M= 3.3
s= 0.0
-----------------------
k= 2
M= 4.15
s= 1.4449999999999996
-----------------------
k= 3
M= 5.166666666666667
s= 7.646666666666666
-----------------------
k= 4
M= 6.875
s= 42.6675
-----------------------
k= 5
M= 6.3
s= 49.279999999999994
-----------------------
k= 6
M= 6.25
s= 49.355
-----------------------
(6.828571428571428, 10.56904761904762)
Calculons et comparons en fait la moyenne et la variance avec numpy. L'argument mot-clé ddof = 1 sert à calculer comme une distribution non biaisée.
import numpy as np
x_arr = np.array(x)
print((x_arr.mean(), x_arr.var(ddof=1)))
Résultat de sortie
(6.828571428571428, 10.569047619047621)
J'ai pu calculer séquentiellement comme ça. En fait, je ne pense pas que python récursif soit plus rapide que numpy, donc je ne pense pas avoir une chance de l'implémenter moi-même, mais bien qu'il soit excellent comme algorithme de calcul de distribution, il n'est pas bien connu, alors laissez-moi l'écrire sous forme d'article. Je l'ai reçu. J'espère que cela sera utile à tout le monde.
Accurately computing running variance https://www.johndcook.com/blog/standard_deviation/
Recommended Posts