Algorithmen zur Berechnung der Varianz - Algorithms for calculating variance

Algorithmen zur Berechnung der Varianz spielen in der Computerstatistik eine große Rolle . Eine Hauptschwierigkeit beim Entwurf guter Algorithmen für dieses Problem besteht darin, dass Formeln für die Varianz Quadratsummen enthalten können, was zu numerischer Instabilität sowie zu arithmetischem Überlauf bei großen Werten führen kann.

Naiver Algorithmus

Eine Formel zur Berechnung der Varianz einer Gesamtpopulation der Größe N lautet:

${\displaystyle \sigma^{2}={\overline {(x^{2})}}-{\bar {x}}^{2}=\displaystyle {\frac {\sum _{i=1} ^{N}x_{i}^{2}-(\sum_{i=1}^{N}x_{i})^{2}/N}{N}}.\!}$

Unter Verwendung der Bessel-Korrektur zur Berechnung einer unverzerrten Schätzung der Populationsvarianz aus einer endlichen Stichprobe von n Beobachtungen lautet die Formel:

${\displaystyle s^{2}=\left({\frac {\sum_{i=1}^{n}x_{i}^{2}}{n}}-\left({\frac {\ Summe _{i=1}^{n}x_{i}}{n}}\right)^{2}\right)\cdot {\frac {n}{n-1}}.\!}$

Daher ist ein naiver Algorithmus zur Berechnung der geschätzten Varianz wie folgt gegeben:

Dieser Algorithmus kann leicht angepasst werden, um die Varianz einer endlichen Population zu berechnen: einfach durch N statt n  − 1 in der letzten Zeile dividieren .

Da SumSq und (Sum×Sum)/ n sehr ähnliche Zahlen sein können, kann die Aufhebung dazu führen, dass die Genauigkeit des Ergebnisses viel geringer ist als die inhärente Genauigkeit der zur Durchführung der Berechnung verwendeten Gleitkommaarithmetik . Daher sollte dieser Algorithmus in der Praxis nicht verwendet werden, und es wurden mehrere alternative, numerisch stabile Algorithmen vorgeschlagen. Dies ist besonders schlimm, wenn die Standardabweichung relativ zum Mittelwert klein ist.

Berechnen von verschobenen Daten

Die Varianz ist invariant in Bezug auf Änderungen eines Standortparameters , eine Eigenschaft, die verwendet werden kann, um die katastrophale Auslöschung in dieser Formel zu vermeiden.

${\displaystyle \operatorname {Var} (XK)=\operatorname {Var} (X).}$

mit beliebiger Konstante, was zu der neuen Formel führt ${\displaystyle K}$

${\displaystyle \sigma^{2}=\displaystyle {\frac {\sum_{i=1}^{n}(x_{i}-K)^{2}-(\sum_{i=1} ^{n}(x_{i}-K))^{2}/n}{n-1}}.\!}$

je näher der Mittelwert ist, desto genauer ist das Ergebnis, aber nur die Auswahl eines Wertes innerhalb des Probenbereichs garantiert die gewünschte Stabilität. Wenn die Werte klein sind, gibt es keine Probleme mit der Summe ihrer Quadrate, im Gegenteil, wenn sie groß sind, bedeutet dies zwangsläufig, dass auch die Varianz groß ist. In jedem Fall ist der zweite Term in der Formel immer kleiner als der erste, daher kann keine Aufhebung erfolgen. ${\displaystyle K}$${\displaystyle (x_{i}-K)}$

Wenn nur das erste Beispiel genommen wird, kann der Algorithmus in der Programmiersprache Python geschrieben werden als ${\displaystyle K}$

def shifted_data_variance(data):
    if len(data) < 2:
        return 0.0
    K = data[0]
    n = Ex = Ex2 = 0.0
    for x in data:
        n = n + 1
        Ex += x - K
        Ex2 += (x - K) * (x - K)
    variance = (Ex2 - (Ex * Ex) / n) / (n - 1)
    # use n instead of (n-1) if want to compute the exact variance of the given data
    # use (n-1) if data are samples of a larger population
    return variance

Diese Formel erleichtert auch die inkrementelle Berechnung, die ausgedrückt werden kann als

K = n = Ex = Ex2 = 0.0

def add_variable(x):
    global K, n, Ex, Ex2
    if n == 0:
        K = x
    n += 1
    Ex += x - K
    Ex2 += (x - K) * (x - K)

def remove_variable(x):
    global K, n, Ex, Ex2
    n -= 1
    Ex -= x - K
    Ex2 -= (x - K) * (x - K)

def get_mean():
    global K, n, Ex
    return K + Ex / n

def get_variance():
    global n, Ex, Ex2
    return (Ex2 - (Ex * Ex) / n) / (n - 1)

Zwei-Pass-Algorithmus

Ein alternativer Ansatz, der eine andere Formel für die Varianz verwendet, berechnet zuerst den Stichprobenmittelwert,

${\displaystyle {\bar{x}}={\frac {\sum_{j=1}^{n}x_{j}}{n}},}$

und berechnet dann die Summe der Quadrate der Differenzen vom Mittelwert,

${\displaystyle {\text{Stichprobenvarianz}}=s^{2}=\displaystyle {\frac {\sum _{i=1}^{n}(x_{i}-{\bar {x}}) ^{2}}{n-1}},\!}$

wobei s die Standardabweichung ist. Dies ist durch den folgenden Code gegeben:

def two_pass_variance(data):
    n = sum1 = sum2 = 0

    for x in data:
        n += 1
        sum1 += x

    mean = sum1 / n

    for x in data:
        sum2 += (x - mean) * (x - mean)

    variance = sum2 / (n - 1)
    return variance

Dieser Algorithmus ist numerisch stabil, wenn n klein ist. Die Ergebnisse dieser beiden einfachen Algorithmen ("naiv" und "two-pass") können jedoch übermäßig von der Reihenfolge der Daten abhängen und können bei sehr großen Datensätzen aufgrund wiederholter Rundungsfehler bei der Akkumulation der Summen. Techniken wie kompensierte Summation können verwendet werden, um diesen Fehler bis zu einem gewissen Grad zu bekämpfen.

Welfords Online-Algorithmus

Es ist oft nützlich, die Varianz in einem einzigen Durchgang berechnen zu können, indem jeder Wert nur einmal überprüft wird; zum Beispiel, wenn die Daten ohne genügend Speicherplatz gesammelt werden, um alle Werte zu behalten, oder wenn die Kosten für den Speicherzugriff die der Berechnung überwiegen. Für einen solchen Online-Algorithmus ist eine Rekursionsbeziehung zwischen Größen erforderlich, aus der die erforderliche Statistik numerisch stabil berechnet werden kann. ${\displaystyle x_{i}}$

Die folgenden Formeln können verwendet werden, um den Mittelwert und die (geschätzte) Varianz der Sequenz für ein zusätzliches Element x _n zu aktualisieren . Hier bezeichnet die Probe Mittelwert der ersten n Proben , deren vorgespannte Stichprobenvarianz , und ihre unvoreingenommenen Stichprobenvarianz . ${\displaystyle {\overline {x}}_{n}={\frac {1}{n}}\sum _{i=1}^{n}x_{i}}$${\displaystyle (x_{1},\dots,x_{n})}$${\displaystyle \sigma _{n}^{2}={\frac {1}{n}}\sum _{i=1}^{n}\left(x_{i}-{\overline {x}) }_{n}\right)^{2}}$${\displaystyle s_{n}^{2}={\frac {1}{n-1}}\sum _{i=1}^{n}\left(x_{i}-{\overline {x}) }_{n}\right)^{2}}$

${\displaystyle {\bar {x}}_{n}={\frac {(n-1)\,{\bar {x}}_{n-1}+x_{n}}{n}}= {\bar {x}}_{n-1}+{\frac {x_{n}-{\bar {x}}_{n-1}}{n}}\!}$

${\displaystyle \sigma_{n}^{2}={\frac {(n-1)\,\sigma_{n-1}^{2}+(x_{n}-{\bar {x}) }_{n-1})(x_{n}-{\bar{x}}_{n})}{n}}=\sigma_{n-1}^{2}+{\frac {( x_{n}-{\bar{x}}_{n-1})(x_{n}-{\bar{x}}_{n})-\sigma_{n-1}^{2} }{n}}.}$

${\displaystyle s_{n}^{2}={\frac {n-2}{n-1}}\,s_{n-1}^{2}+{\frac {(x_{n}-{ \bar{x}}_{n-1})^{2}}{n}}=s_{n-1}^{2}+{\frac {(x_{n}-{\bar {x}) }_{n-1})^{2}}{n}}-{\frac {s_{n-1}^{2}}{n-1}},\quad n>1}$

Diese Formeln leiden unter numerischer Instabilität, da sie wiederholt eine kleine Zahl von einer großen Zahl subtrahieren, die mit n skaliert . Eine bessere Größe für die Aktualisierung ist die Summe der Quadrate der Differenzen vom aktuellen Mittelwert , hier bezeichnet : ${\displaystyle \textstyle \sum_{i=1}^{n}(x_{i}-{\bar {x}}_{n})^{2}}$${\displaystyle M_{2,n}}$

${\displaystyle {\begin{ausgerichtet}M_{2,n}&=M_{2,n-1}+(x_{n}-{\bar {x}}_{n-1})(x_{n }-{\bar {x}}_{n})\\[4pt]\sigma_{n}^{2}&={\frac {M_{2,n}}{n}}\\[4pt ]s_{n}^{2}&={\frac {M_{2,n}}{n-1}}\end{ausgerichtet}}}$

Dieser Algorithmus wurde von Welford gefunden und gründlich analysiert. Es ist auch üblich, und zu bezeichnen . ${\displaystyle M_{k}={\bar {x}}_{k}}$${\displaystyle S_{k}=M_{2,k}}$

Eine beispielhafte Python-Implementierung für den Welford-Algorithmus ist unten angegeben.

# For a new value newValue, compute the new count, new mean, the new M2.
# mean accumulates the mean of the entire dataset
# M2 aggregates the squared distance from the mean
# count aggregates the number of samples seen so far
def update(existingAggregate, newValue):
    (count, mean, M2) = existingAggregate
    count += 1
    delta = newValue - mean
    mean += delta / count
    delta2 = newValue - mean
    M2 += delta * delta2
    return (count, mean, M2)

# Retrieve the mean, variance and sample variance from an aggregate
def finalize(existingAggregate):
    (count, mean, M2) = existingAggregate
    if count < 2:
        return float("nan")
    else:
        (mean, variance, sampleVariance) = (mean, M2 / count, M2 / (count - 1))
        return (mean, variance, sampleVariance)

Dieser Algorithmus ist viel weniger anfällig für Präzisionsverluste aufgrund einer katastrophalen Aufhebung , ist jedoch möglicherweise aufgrund der Divisionsoperation innerhalb der Schleife nicht so effizient. Für einen besonders robusten Zweidurchgangsalgorithmus zur Berechnung der Varianz kann man zuerst eine Schätzung des Mittelwerts berechnen und subtrahieren und dann diesen Algorithmus auf die Residuen anwenden.

Der folgende parallele Algorithmus veranschaulicht, wie mehrere online berechnete Statistiken zusammengeführt werden.

Gewichteter inkrementeller Algorithmus

Der Algorithmus kann erweitert werden, um ungleiche Stichprobengewichte zu handhaben, indem der einfache Zähler n durch die Summe der bisher gesehenen Gewichte ersetzt wird. West (1979) schlägt diesen inkrementellen Algorithmus vor :

def weighted_incremental_variance(data_weight_pairs):
    w_sum = w_sum2 = mean = S = 0

    for x, w in data_weight_pairs:
        w_sum = w_sum + w
        w_sum2 = w_sum2 + w * w
        mean_old = mean
        mean = mean_old + (w / w_sum) * (x - mean_old)
        S = S + w * (x - mean_old) * (x - mean)

    population_variance = S / w_sum
    # Bessel's correction for weighted samples
    # Frequency weights
    sample_frequency_variance = S / (w_sum - 1)
    # Reliability weights
    sample_reliability_variance = S / (w_sum - w_sum2 / w_sum)

Paralleler Algorithmus

Chanet al. Beachten Sie, dass der oben beschriebene Online-Algorithmus von Welford ein Sonderfall eines Algorithmus ist, der zum Kombinieren beliebiger Mengen und funktioniert : ${\displaystyle A}$${\displaystyle B}$

${\displaystyle {\begin{ausgerichtet}n_{AB}&=n_{A}+n_{B}\\\delta &={\bar {x}}_{B}-{\bar {x}}_ {A}\\{\bar{x}}_{AB}&={\bar{x}}_{A}+\delta \cdot {\frac {n_{B}}{n_{AB}}} \\M_{2,AB}&=M_{2,A}+M_{2,B}+\delta ^{2}\cdot {\frac{n_{A}n_{B}}{n_{AB} }}\\Ende{ausgerichtet}}}$.

Dies kann nützlich sein, wenn beispielsweise diskreten Teilen der Eingabe mehrere Verarbeitungseinheiten zugewiesen werden können.

Chans Methode zur Schätzung des Mittelwerts ist numerisch instabil, wenn und beide groß sind, da der numerische Fehler in nicht so verkleinert wird, wie es im Fall der Fall ist. In solchen Fällen bevorzugen Sie . ${\displaystyle n_{A}\approx n_{B}}$${\displaystyle \delta ={\bar{x}}_{B}-{\bar{x}}_{A}}$${\displaystyle n_{B}=1}$${\displaystyle {\bar {x}}_{AB}={\frac {n_{A}{\bar {x}}_{A}+n_{B}{\bar {x}}_{B} }{schnappen}}}}$

def parallel_variance(n_a, avg_a, M2_a, n_b, avg_b, M2_b):
    n = n_a + n_b
    delta = avg_b - avg_a
    M2 = M2_a + M2_b + delta ** 2 * n_a * n_b / n
    var_ab = M2 / (n - 1)
    return var_ab

Dies kann verallgemeinert werden, um Parallelisierung mit AVX , mit GPUs und Computerclustern und Kovarianz zu ermöglichen.

Beispiel

Angenommen, alle Gleitkommaoperationen verwenden die Standard- IEEE 754- Arithmetik mit doppelter Genauigkeit . Betrachten Sie die Stichprobe (4, 7, 13, 16) aus einer unendlichen Grundgesamtheit. Basierend auf dieser Stichprobe beträgt der geschätzte Populationsmittelwert 10 und die unverzerrte Schätzung der Populationsvarianz 30. Sowohl der naive Algorithmus als auch der Two-Pass-Algorithmus berechnen diese Werte korrekt.

Betrachten Sie als nächstes die Stichprobe ( 10 ⁸  + 4 , 10 ⁸  + 7 , 10 ⁸  + 13 , 10 ⁸  + 16 ), die dieselbe geschätzte Varianz wie die erste Stichprobe ergibt. Der Algorithmus mit zwei Durchläufen berechnet diese Varianzschätzung korrekt, aber der naive Algorithmus gibt 29.33333333333332 statt 30 zurück.

Während dieser Genauigkeitsverlust tolerierbar und als kleiner Fehler des naiven Algorithmus angesehen werden kann, macht eine weitere Erhöhung des Offsets den Fehler katastrophal. Betrachten Sie die Stichprobe ( 10 ⁹  + 4 , 10 ⁹  + 7 , 10 ⁹  + 13 , 10 ⁹  + 16 ). Auch hier wird die geschätzte Populationsvarianz von 30 vom Two-Pass-Algorithmus korrekt berechnet, aber der naive Algorithmus berechnet sie jetzt als −170,66666666666666. Dies ist ein ernstes Problem bei naiven Algorithmen und ist auf eine katastrophale Aufhebung bei der Subtraktion zweier ähnlicher Zahlen in der Endstufe des Algorithmus zurückzuführen.

Statistik höherer Ordnung

Terriberry erweitert Chans Formeln auf die Berechnung des dritten und vierten zentralen Moments , die zum Beispiel bei der Schätzung von Schiefe und Kurtosis benötigt werden :

${\displaystyle {\begin{ausgerichtet}M_{3,X}=M_{3,A}+M_{3,B}&{}+\delta ^{3}{\frac {n_{A}n_{B }(n_{A}-n_{B})}{n_{X}^{2}}}+3\delta {\frac {n_{A}M_{2,B}-n_{B}M_{2 ,A}}{n_{X}}}\\[6pt]M_{4,X}=M_{4,A}+M_{4,B}&{}+\delta ^{4}{\frac { n_{A}n_{B}\left(n_{A}^{2}-n_{A}n_{B}+n_{B}^{2}\right)}{n_{X}^{3} }}\\[6pt]&{}+6\delta ^{2}{\frac {n_{A}^{2}M_{2,B}+n_{B}^{2}M_{2,A }}{n_{X}^{2}}}+4\delta {\frac{n_{A}M_{3,B}-n_{B}M_{3,A}}{n_{X}}} \end{ausgerichtet}}}$

Hier sind wieder die Summen der Potenzen der Differenzen vom Mittelwert , die ergeben ${\displaystyle M_{k}}$${\displaystyle \sum (x-{\overline {x}})^{k}}$

${\displaystyle {\begin{ausgerichtet}&{\text{Schiefe}}=g_{1}={\frac {{\sqrt {n}}M_{3}}{M_{2}^{3/2} }},\\[4pt]&{\text{kurtosis}}=g_{2}={\frac {nM_{4}}{M_{2}^{2}}}-3.\end{aligned} }}$

Für den inkrementellen Fall (dh ) vereinfacht sich dies zu: ${\displaystyle B=\{x\}}$

${\displaystyle {\begin{aligned}\delta &=xm\\[5pt]m'&=m+{\frac {\delta }{n}}\\[5pt]M_{2}'&=M_{2 }+\delta ^{2}{\frac {n-1}{n}}\\[5pt]M_{3}'&=M_{3}+\delta ^{3}{\frac {(n- 1)(n-2)}{n^{2}}}-{\frac {3\delta M_{2}}{n}}\\[5pt]M_{4}'&=M_{4}+ {\frac {\delta^{4}(n-1)(n^{2}-3n+3)}{n^{3}}}+{\frac {6\delta^{2}M_{2 }}{n^{2}}}-{\frac {4\delta M_{3}}{n}}\end{ausgerichtet}}}$

Durch Beibehalten des Wertes wird nur eine Divisionsoperation benötigt und die Statistik höherer Ordnung kann somit mit geringen zusätzlichen Kosten berechnet werden. ${\displaystyle \delta /n}$

Ein Beispiel für den wie beschrieben implementierten Online-Algorithmus für Kurtosis ist:

def online_kurtosis(data):
    n = mean = M2 = M3 = M4 = 0

    for x in data:
        n1 = n
        n = n + 1
        delta = x - mean
        delta_n = delta / n
        delta_n2 = delta_n * delta_n
        term1 = delta * delta_n * n1
        mean = mean + delta_n
        M4 = M4 + term1 * delta_n2 * (n*n - 3*n + 3) + 6 * delta_n2 * M2 - 4 * delta_n * M3
        M3 = M3 + term1 * delta_n * (n - 2) - 3 * delta_n * M2
        M2 = M2 + term1

    # Note, you may also calculate variance using M2, and skewness using M3
    # Caution: If all the inputs are the same, M2 will be 0, resulting in a division by 0.
    kurtosis = (n * M4) / (M2 * M2) - 3
    return kurtosis

Pébaÿ erweitert diese Ergebnisse weiter auf Zentralmomente beliebiger Ordnung für den inkrementellen und den paarweisen Fall, und anschließend Pébaÿ et al. für gewichtete und zusammengesetzte Momente. Dort findet man auch ähnliche Formeln für die Kovarianz .

Choi und Sweetman bieten zwei alternative Methoden zum Berechnen der Schiefe und Kurtosis an, von denen jede in bestimmten Anwendungen erheblichen Computerspeicherbedarf und CPU-Zeit einsparen kann. Der erste Ansatz besteht darin, die statistischen Momente zu berechnen, indem die Daten in Bins aufgeteilt werden und dann die Momente aus der Geometrie des resultierenden Histogramms berechnet werden , was effektiv zu einem Algorithmus mit einem Durchgang für höhere Momente wird. Ein Vorteil besteht darin, dass die Berechnungen des statistischen Moments mit beliebiger Genauigkeit durchgeführt werden können, so dass die Berechnungen auf die Genauigkeit beispielsweise des Datenspeicherformats oder der ursprünglichen Messhardware abgestimmt werden können. Ein relatives Histogramm einer Zufallsvariable kann auf herkömmliche Weise konstruiert werden: Der Bereich potenzieller Werte wird in Bins unterteilt und die Anzahl der Vorkommen innerhalb jeder Bin wird gezählt und so aufgetragen, dass die Fläche jedes Rechtecks ​​dem Anteil der Sample-Werte entspricht in diesem Behälter:

${\displaystyle H(x_{k})={\frac {h(x_{k})}{A}}}$

wobei und die Häufigkeit und die relative Häufigkeit bei bin darstellen und die Gesamtfläche des Histogramms ist. Nach dieser Normierung können aus dem relativen Histogramm die Rohmomente und Zentralmomente von berechnet werden: ${\displaystyle h(x_{k})}$${\displaystyle H(x_{k})}$${\displaystyle x_{k}}$${\displaystyle A=\sum_{k=1}^{K}h(x_{k})\,\Delta x_{k}}$${\displaystyle n}$${\displaystyle x(t)}$

${\displaystyle m_{n}^{(h)}=\sum_{k=1}^{K}x_{k}^{n}H(x_{k})\,\Updelta x_{k}= {\frac {1}{A}}\sum_{k=1}^{K}x_{k}^{n}h(x_{k})\,\Updelta x_{k}}$

${\displaystyle \theta_{n}^{(h)}=\sum_{k=1}^{K}{\Groß (}x_{k}-m_{1}^{(h)}{\ Groß)}^{n}\,H(x_{k})\,\Delta x_{k}={\frac {1}{A}}\sum _{k=1}^{K}{\Groß (}x_{k}-m_{1}^{(h)}{\Groß)}^{n}h(x_{k})\,\Delta x_{k}}$

wobei der hochgestellte Index angibt, dass die Momente aus dem Histogramm berechnet werden. Für eine konstante Bin-Breite können diese beiden Ausdrücke mit vereinfacht werden : ${\displaystyle ^{(h)}}$${\displaystyle \Updelta x_{k}=\Updelta x}$${\displaystyle I=A/\Delta x}$

${\displaystyle m_{n}^{(h)}={\frac {1}{I}}\sum _{k=1}^{K}x_{k}^{n}\,h(x_{ k})}$

${\displaystyle \theta_{n}^{(h)}={\frac {1}{I}}\sum _{k=1}^{K}{\Big (}x_{k}-m_{ 1}^{(h)}{\Groß)}^{n}h(x_{k})}$

Der zweite Ansatz von Choi und Sweetman ist eine analytische Methodik, um statistische Momente aus einzelnen Segmenten einer Zeitgeschichte so zu kombinieren, dass die resultierenden Gesamtmomente die der gesamten Zeitgeschichte sind. Diese Methodik könnte für die parallele Berechnung statistischer Momente mit anschließender Kombination dieser Momente oder für die Kombination von statistischen Momenten, die zu sequentiellen Zeiten berechnet werden, verwendet werden.

Wenn Sätze von statistischen Momenten bekannt sind: für , dann kann jeder durch die äquivalenten Rohmomente ausgedrückt werden : ${\displaystyle Q}$${\displaystyle (\gamma_{0,q},\mu_{q},\sigma_{q}^{2},\alpha_{3,q},\alpha_{4,q})\ Quad }$${\displaystyle q=1,2,\ldots,Q}$${\displaystyle \gamma_{n}}$${\displaystyle n}$

${\displaystyle \gamma_{n,q}=m_{n,q}\gamma_{0,q}\qquad\quad {\textrm {für}}\quad n=1,2,3,4\quad {\text{ und }}\quad q=1,2,\dots,Q}$

wobei im Allgemeinen die Dauer des Zeitverlaufs oder die Anzahl der Punkte, wenn konstant ist, angenommen wird. ${\displaystyle \gamma_{0,q}}$${\displaystyle q^{th}}$${\displaystyle \Delta t}$

Der Vorteil, die statistischen Momente in Form von auszudrücken, besteht darin, dass die Sätze durch Addition kombiniert werden können und es keine Obergrenze für den Wert von gibt . ${\displaystyle\gamma}$${\displaystyle Q}$${\displaystyle Q}$

${\displaystyle \gamma_{n,c}=\sum_{q=1}^{Q}\gamma_{n,q}\quad\quad {\text{für}}n=0,1,2 ,3,4}$

wobei der Index den verketteten Zeitverlauf oder kombinierte . Diese kombinierten Werte von können dann umgekehrt in Rohmomente umgewandelt werden, die die vollständige verkettete Zeitgeschichte darstellen ${\displaystyle_{c}}$${\displaystyle\gamma}$${\displaystyle\gamma}$

${\displaystyle m_{n,c}={\frac {\gamma_{n,c}}{\gamma_{0,c}}}\quad {\text{für}}n=1,2,3 ,4}$

Bekannte Beziehungen zwischen den Rohmomenten ( ) und den Zentralmomenten ( ) werden dann verwendet, um die Zentralmomente des verketteten Zeitverlaufs zu berechnen. Aus den zentralen Momenten werden schließlich die statistischen Momente der verketteten Geschichte berechnet: ${\displaystyle m_{n}}$${\displaystyle \theta_{n}=\operatorname {E} [(x-\mu)^{n}])}$

${\displaystyle \mu_{c}=m_{1,c}\qquad\sigma_{c}^{2}=\theta_{2,c}\qquad\alpha_{3,c}={\ frac {\theta_{3,c}}{\sigma_{c}^{3}}}\qquad\alpha_{4,c}={\frac{\theta_{4,c}}{\ Sigma _{c}^{4}}}-3}$

Kovarianz

Sehr ähnliche Algorithmen können verwendet werden, um die Kovarianz zu berechnen .

Naiver Algorithmus

Der naive Algorithmus lautet:

${\displaystyle \operatorname {Cov} (X,Y)=\displaystyle {\frac {\sum _{i=1}^{n}x_{i}y_{i}-(\sum _{i=1} ^{n}x_{i})(\sum_{i=1}^{n}y_{i})/n}{n}}.\!}$

Für den obigen Algorithmus könnte man den folgenden Python-Code verwenden:

def naive_covariance(data1, data2):
    n = len(data1)
    sum12 = 0
    sum1 = sum(data1)
    sum2 = sum(data2)

    for i1, i2 in zip(data1, data2):
        sum12 += i1 * i2

    covariance = (sum12 - sum1 * sum2 / n) / n
    return covariance

Mit Schätzung des Mittelwerts

Was die Varianz betrifft, so ist auch die Kovarianz zweier Zufallsvariablen verschiebungsinvariant, also gegeben zwei konstante Werte und kann geschrieben werden: ${\displaystyle k_{x}}$${\displaystyle k_{y},}$

${\displaystyle \operatorname {Cov} (X,Y)=\operatorname {Cov} (X-k_{x},Y-k_{y})=\displaystyle {\frac {\sum _{i=1}^ {n}(x_{i}-k_{x})(y_{i}-k_{y})-(\sum_{i=1}^{n}(x_{i}-k_{x}) )(\sum_{i=1}^{n}(y_{i}-k_{y}))/n}{n}}.\!}$

und wieder wird die Wahl eines Wertes innerhalb des Wertebereichs die Formel gegen katastrophale Annullierung stabilisieren und sie robuster gegen große Summen machen. Nimmt man den ersten Wert jedes Datensatzes, kann der Algorithmus wie folgt geschrieben werden:

def shifted_data_covariance(data_x, data_y):
    n = len(data_x)
    if n < 2:
        return 0
    kx = data_x[0]
    ky = data_y[0]
    Ex = Ey = Exy = 0
    for ix, iy in zip(data_x, data_y):
        Ex += ix - kx
        Ey += iy - ky
        Exy += (ix - kx) * (iy - ky)
    return (Exy - Ex * Ey / n) / n

Zwei-Pass

Der Two-Pass-Algorithmus berechnet zuerst die Stichprobenmittelwerte und dann die Kovarianz:

${\displaystyle {\bar {x}}=\displaystyle \sum _{i=1}^{n}x_{i}/n}$

${\displaystyle {\bar {y}}=\displaystyle \sum _{i=1}^{n}y_{i}/n}$

${\displaystyle \operatorname {Cov} (X,Y)=\displaystyle {\frac {\sum _{i=1}^{n}(x_{i}-{\bar {x}})(y_{i }-{\bar{y}})}{n}}.\!}$

Der Algorithmus mit zwei Durchläufen kann wie folgt geschrieben werden:

def two_pass_covariance(data1, data2):
    n = len(data1)

    mean1 = sum(data1) / n
    mean2 = sum(data2) / n

    covariance = 0

    for i1, i2 in zip(data1, data2):
        a = i1 - mean1
        b = i2 - mean2
        covariance += a * b / n
    return covariance

Eine etwas genauere kompensierte Version führt den vollständigen naiven Algorithmus an den Residuen durch. Die Endsummen und sollten null sein, aber der zweite Durchgang kompensiert jeden kleinen Fehler. ${\displaystyle \textstyle \sum x_{i}}$${\displaystyle \textstyle \sum y_{i}}$

Online

Es existiert ein stabiler One-Pass-Algorithmus, ähnlich dem Online-Algorithmus zur Berechnung der Varianz, der das Co-Moment berechnet : ${\displaystyle \textstyle C_{n}=\sum _{i=1}^{n}(x_{i}-{\bar {x}}_{n})(y_{i}-{\bar { j}}_{n})}$

${\displaystyle {\begin{alignedat}{2}{\bar {x}}_{n}&={\bar {x}}_{n-1}&\,+\,&{\frac {x_ {n}-{\bar {x}}_{n-1}}{n}}\\[5pt]{\bar {y}}_{n}&={\bar {y}}_{n -1}&\,+\,&{\frac {y_{n}-{\bar {y}}_{n-1}}{n}}\\[5pt]C_{n}&=C_{ n-1}&\,+\,&(x_{n}-{\bar{x}}_{n})(y_{n}-{\bar{y}}_{n-1})\ \[5pt]&=C_{n-1}&\,+\,&(x_{n}-{\bar {x}}_{n-1})(y_{n}-{\bar {y }}_{n})\end{alignedat}}}$

Die scheinbare Asymmetrie in dieser letzten Gleichung ist darauf zurückzuführen, dass also beide Aktualisierungsterme gleich sind . Eine noch größere Genauigkeit kann erreicht werden, indem zuerst die Mittelwerte berechnet werden und dann der stabile Einmal-Algorithmus für die Residuen verwendet wird. ${\displaystyle \textstyle (x_{n}-{\bar {x}}_{n})={\frac {n-1}{n}}(x_{n}-{\bar {x}}_ {n-1})}$${\displaystyle \textstyle {\frac {n-1}{n}}(x_{n}-{\bar {x}}_{n-1})(y_{n}-{\bar {y}} _{n-1})}$

Damit kann die Kovarianz berechnet werden als

${\displaystyle {\begin{ausgerichtet}\operatorname {Cov} _{N}(X,Y)={\frac {C_{N}}{N}}&={\frac {\operatorname {Cov} _{ N-1}(X,Y)\cdot (N-1)+(x_{n}-{\bar{x}}_{n})(y_{n}-{\bar{y}}_{ n-1})}{N}}\\&={\frac {\operatorname {Cov} _{N-1}(X,Y)\cdot (N-1)+(x_{n}-{\ bar {x}}_{n-1})(y_{n}-{\bar {y}}_{n})}{N}}\\&={\frac {\operatorname {Cov} _{ N-1}(X,Y)\cdot (N-1)+{\frac {N-1}{N}}(x_{n}-{\bar{x}}_{n-1})( y_{n}-{\bar {y}}_{n-1})}{N}}\\&={\frac {\operatorname {Cov} _{N-1}(X,Y)\cdot (N-1)+{\frac{N}{N-1}}(x_{n}-{\bar{x}}_{n})(y_{n}-{\bar{y}}_ {n})}{N}}.\end{aligned}}}$

def online_covariance(data1, data2):
    meanx = meany = C = n = 0
    for x, y in zip(data1, data2):
        n += 1
        dx = x - meanx
        meanx += dx / n
        meany += (y - meany) / n
        C += dx * (y - meany)

    population_covar = C / n
    # Bessel's correction for sample variance
    sample_covar = C / (n - 1)

Eine kleine Modifikation kann auch vorgenommen werden, um die gewichtete Kovarianz zu berechnen:

def online_weighted_covariance(data1, data2, data3):
    meanx = meany = 0
    wsum = wsum2 = 0
    C = 0
    for x, y, w in zip(data1, data2, data3):
        wsum += w
        wsum2 += w * w
        dx = x - meanx
        meanx += (w / wsum) * dx
        meany += (w / wsum) * (y - meany)
        C += w * dx * (y - meany)

    population_covar = C / wsum
    # Bessel's correction for sample variance
    # Frequency weights
    sample_frequency_covar = C / (wsum - 1)
    # Reliability weights
    sample_reliability_covar = C / (wsum - wsum2 / wsum)

Ebenso gibt es eine Formel zum Kombinieren der Kovarianzen zweier Mengen, die verwendet werden kann, um die Berechnung zu parallelisieren:

${\displaystyle C_{X}=C_{A}+C_{B}+({\bar {x}}_{A}-{\bar {x}}_{B})({\bar {y} }_{A}-{\bar{y}}_{B})\cdot {\frac {n_{A}n_{B}}{n_{X}}}.}$

Gewichtete Batch-Version

Es gibt auch eine Version des gewichteten Online-Algorithmus, die Batch-Updates durchführt: Bezeichne die Gewichte und schreibe ${\displaystyle w_{1},\dots w_{N}}$

${\displaystyle {\begin{alignedat}{2}{\bar {x}}_{n+k}&={\bar {x}}_{n}&\,+\,&{\frac {\ Summe _{i=n+1}^{n+k}w_{i}(x_{i}-{\bar {x}}_{n})}{\sum _{i=1}^{n +k}w_{i}}}\\{\bar{y}}_{n+k}&={\bar {y}}_{n}&\,+\,&{\frac {\sum _{i=n+1}^{n+k}w_{i}(y_{i}-{\bar {y}}_{n})}{\sum _{i=1}^{n+ k}w_{i}}}\\C_{n+k}&=C_{n}&\,+\,&\sum _{i=n+1}^{n+k}w_{i}( x_{i}-{\bar {x}}_{n+k})(y_{i}-{\bar {y}}_{n})\\&=C_{n}&\,+\ ,&\sum_{i=n+1}^{n+k}w_{i}(x_{i}-{\bar {x}}_{n})(y_{i}-{\bar { y}}_{n+k})\\\end{alignedat}}}$

Die Kovarianz kann dann berechnet werden als

${\displaystyle \operatorname {Cov} _{N}(X,Y)={\frac {C_{N}}{\sum _{i=1}^{N}w_{i}}}}$

Siehe auch

• Kahan-Summierungsalgorithmus
• Quadratische Abweichungen vom Mittelwert
• Yamartino-Methode

Verweise

Externe Links

• Weisstein, Eric W. "Sample Variance Computation" . MathWorld .