Autokorrelation zwischen Python und Julia

19

Ich versuche, mit Julia eine Autokorrelation durchzuführen und sie mit Pythons Ergebnis zu vergleichen. Wie kommt es, dass sie unterschiedliche Ergebnisse liefern?

Julia Code

using StatsBase

t = range(0, stop=10, length=10)
test_data = sin.(exp.(t.^2))

acf = StatsBase.autocor(test_data)

gibt

10-element Array{Float64,1}:
  1.0                   
  0.13254954979179642   
 -0.2030283419321465    
  0.00029587850872956104
 -0.06629381497277881   
  0.031309038331589614  
 -0.16633393452504994   
 -0.08482388975165675   
  0.0006905628640697538 
 -0.1443650483145533

Python-Code

from statsmodels.tsa.stattools import acf
import numpy as np

t = np.linspace(0,10,10)
test_data = np.sin(np.exp(t**2))

acf_result = acf(test_data)

gibt

array([ 1.        ,  0.14589844, -0.10412699,  0.07817509, -0.12916543,
       -0.03469143, -0.129255  , -0.15982435, -0.02067688, -0.14633346])
Ross Mariano
quelle
1
Drucken Sie die Testdaten in beiden Fällen aus
Mad Physicist

Antworten:

26

Das liegt daran, dass Sie test_dataanders sind:

Python:

array([ 0.84147098, -0.29102733,  0.96323736,  0.75441021, -0.37291918,
        0.85600145,  0.89676529, -0.34006519, -0.75811102, -0.99910501])

Julia:

[0.8414709848078965, -0.2910273263243299, 0.963237364649543, 0.7544102058854344,
 -0.3729191776326039, 0.8560014512776061, 0.9841238290665676, 0.1665709194875013,
 -0.7581110212957692, -0.9991050130774393]

Dies geschieht, weil Sie sinenorme Zahlen nehmen. Zum Beispiel ist die letzte Zahl t10 exp(10^2)~ 2,7 * 10 ^ 43. Auf dieser Skala betragen die Gleitkomma-Ungenauigkeiten etwa 3 * 10 ^ 9. Wenn also auch das niedrigstwertige Bit für Python und Julia unterschiedlich ist, ist der sinWert weit entfernt.

Tatsächlich können wir die zugrunde liegenden Binärwerte des anfänglichen Arrays untersuchen t. Zum Beispiel unterscheiden sie sich im drittletzten Wert:

Julia:

julia> reinterpret(Int, range(0, stop=10, length=10)[end-2])
4620443017702830535

Python:

>>> import struct
>>> s = struct.pack('>d', np.linspace(0,10,10)[-3])
>>> struct.unpack('>q', s)[0]
4620443017702830536

Wir können in der Tat sehen, dass sie durch genau eine Maschine epsilon nicht übereinstimmen. Und wenn wir Julia verwenden, nehmen Sie sinden von Python erhaltenen Wert:

julia> sin(exp(reinterpret(Float64, 4620443017702830536)^2))
-0.3400651855865199

Wir erhalten den gleichen Wert, den Python hat.

Jakob Nissen
quelle
9

Nur um die Antwort ein wenig zu erweitern (als Antwort hinzufügen, da sie für einen Kommentar zu lang ist). In Julia haben Sie Folgendes:

julia> t = collect(range(0, stop=10, length=10))
10-element Array{Float64,1}:
  0.0               
  1.1111111111111112
  2.2222222222222223
  3.3333333333333335
  4.444444444444445 
  5.555555555555555 
  6.666666666666667 
  7.777777777777778 
  8.88888888888889  
 10.0               

julia> t .- [10*i / 9 for i in 0:9]
10-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

während in Python:

>>> t = np.linspace(0,10,10)
>>> t - [10*i/9 for i in range(10)]
array([0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,
       0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 8.8817842e-16,
       0.0000000e+00, 0.0000000e+00])

und Sie sehen, dass die 8. Zahl in Python eine ungenaue Annäherung an ist 70/9, während Sie in Julia in diesem Fall die Folge der engsten Annäherungen an die 10*i/9Verwendung erhalten Float64.

Es scheint also, dass der Rest dem folgt, was @Jakob Nissen kommentiert hat, weil sich die ursprünglichen Sequenzen unterscheiden.

Die Dinge sind jedoch nicht so einfach. Da expFunktionen in Julia und Python sich ein wenig in dem unterscheiden, was sie produzieren. Siehe Python:

>>> from math import exp
>>> from mpmath import mp
>>> mp.dps = 1000
>>> float(mp.exp((20/3)**2) - exp((20/3)**2))
-1957.096392544307

während in Julia:

julia> setprecision(1000)
1000

julia> Float64(exp(big((20/3)^2)) - exp((20/3)^2))
2138.903607455693

julia> Float64(exp(big((20/3)^2)) - nextfloat(exp((20/3)^2)))
-1957.096392544307

(Sie können überprüfen, ob (20/3)^2dies Float64sowohl in Julia als auch in Python gleich ist).

In diesem Fall ist expPython also etwas genauer als Julia. Daher wird selbst das Korrigieren t(was durch die Verwendung eines Verständnisses in Python anstelle von einfach ist linspace) nicht dazu führen, dass der ACF gleich ist.

Alles in allem ist die Schlussfolgerung, was @Jakob Nissen für so große Werte kommentierte, dass die Ergebnisse stark von den numerischen Ungenauigkeiten beeinflusst werden.

Bogumił Kamiński
quelle