Časová řada diskrétní veličiny

Jak motivace doby velí, zabýval jsem se různými časovými řadami souvisejícími s pandemií. V Itálii jsou řady nejdelší a dostupné po jednotlivých okresech zde. Vybral jsem si jednu z nich...

In [1]:
import numpy as np
%matplotlib inline
from matplotlib import pyplot as pl
zhash=np.load("../TPX/select_commun.npz",allow_pickle=True)#['arr_0']
zhash=zhash.get('arr_0').reshape(1)[0]
In [37]:
cname='Salerno'
#zhash=zhash.get('arr_0').reshape(1)[0]
pl.plot(zhash[cname],'+')
pl.xlabel('days')
pl.ylabel('positives')
pl.title(cname)
Out[37]:
Text(0.5,1,'Salerno')

Poissonovská veličina se řídí rozdělením $$P(j)=\frac{\lambda(x)^j}{j!}\ \exp(-\lambda(x)),$$

kde $\lambda(x)$ je její očekávaná hodnota (závisí na nějaké další veličině $x$).

Pokud pro sadu bodů $x_i,j_i$ ($i=1-N$) máme model $\lambda(x_i)$, logaritmus věrohodnosti pak vychází

$$\log L = \sum_i \log P (j_i) = \sum_i j_i \log(\lambda(x_i)) - \sum \log(j_i !) - \sum \lambda(x_i)$$

Parametry modelu $\lambda(x)$ pak hledáme derivací této funkce (prostřední člen na $\lambda$ nezávisí, můžeme ho dále vynechat). Předpokládejme kvadratickou závislost $\lambda(x|p_k)=p_0+p_1\ x+p_2\ x^2$

Hledáme teď řešení soustavy $$\sum_i \frac{j_i\ x_i^k}{\lambda(x_i|p_k)} = \sum_i x_i^k$$

Analytické řešení se zdá být obtížné (parametry jsou ve jmenovateli, přesněji $n$ jmenovatelích v sumě vlevo), ne-li nemožné, proto použijeme obecnou numerickou minimalizační metodu.

In [3]:
from scipy import optimize as op
j=zhash[cname]
x=np.r_[:len(j)]
pini=[1,5,3]
lig_ml=lambda p:-(j*np.log(abs(np.polyval(p[::-1],x)))).sum()+np.polyval(p[::-1],x).sum()
par1=op.fmin(lig_ml,pini)
par1

Optimization terminated successfully.
         Current function value: -84309.680558
         Iterations: 174
         Function evaluations: 316
Out[3]:
array([37.28135799, 23.01704175, -0.17879442])

Aproximace normálním rozdělením

Můžeme se pokusit nafitovat tutéž funkci pomocí nejmen. čtverců, což odpovídá předpokládu normálního rozdělení dat. Minimalizujeme zde veličinu $$S=(J-A P)^T\ (J-A P),$$ kde $J$ je vektor měření a modelová matice (design matrix) $$A=\left(\begin{matrix} 1 & x_1 & x_1^2 \\ 1 & x_2 & x_2^2 \\ & \dots & \\ \\ 1 & x_n & x_n^2\end{matrix} \right).\label{desmat} $$ Řešení podmínky pro derivaci podle parametrů $P$ je dáno $\hat P = H^{-1} A^T\ J$, kde $H$ je čtvercová matice $A^T\ A$ (Hessian).

Numericky vychází odhad parametrů takto:

In [4]:
amat=np.array([np.ones_like(x),x,x**2]).T
hess=amat.T@amat
par2=np.linalg.inv(hess)@amat.T@j
par2
Out[4]:
array([20.40191232, 25.47768757, -0.23867128])

Výsledek se výrazně liší od prvního odhadu, problém je ovšem ve výrazné korelaci parametrů (koeficientů polynomu). Inverzní Hessian je úměrný kovarianční matici, korelační z ní určíme normalizací diagonály:

In [13]:
cov=np.linalg.inv(hess)
err=np.sqrt(cov.diagonal())
cor=cov/err[np.newaxis,:]/err[:,np.newaxis]
cor
Out[13]:
array([[ 1.        , -0.742103  ,  0.605679  ],
       [-0.742103  ,  1.        , -0.95669715],
       [ 0.605679  , -0.95669715,  1.        ]])

Je vidět, že korelační koeficient dosahuje až 96% - jde tady o antikorelaci, tedy zvyšování jednoho koeficientu je dobře kompenzováno sni

Ve skutečnosti se od sebe modelové funkce příliš neliší - pod 3 procenta rozsahu dat $J$.

In [5]:
model1=amat@par1
model2=amat@par2
max(abs(model2-model1)), max(j)-min(j)
Out[5]:
(16.87944566963366, 588.0)

Korektní nahrazení poisson. rozdělení pomocí normálně rozdělených proměnných musímě vzít v úvahu proměnný rozptyl: $D(J)=E(J)=\lambda(x)$. V odhadu $P$ se tedy objeví váhová matice W, inverzní ke kovarianční $J$, zde diagonální $W_{ii}=1/D(j_i)$

Pak dosadíme $ H=A^T\ W\ A$ do odhadu $\hat P = H^{-1} A^T\ W\ J$ (je dobré si všimnout, že výsledek nezávisí na normalizaci matice W, můžeme nejistoty měření v jednotlivých bodech znát jen relativně). Odhady i modely se od "poissonovského" výsledku liší jen o necelých 5, tedy méně než jedno procento.

In [6]:
wmat=np.eye(len(j))
wmat[x,x]=1/j #diagonala
hess3=amat.T@wmat@amat
par3=np.linalg.inv(hess3)@amat.T@wmat@j
par3
Out[6]:
array([32.37798359, 23.48434807, -0.18862014])
In [8]:
model1=amat@par1
model3=amat@par3
pl.plot(x,model1, label="Poisson")
pl.plot(x,model3, label="LSM")
pl.plot(x,j,'+',label="data")
pl.legend()
max(abs(model3-model1))
Out[8]:
4.903374404943499

Rozdíl mezi oběma fity už na grafu není viditelný.

Posun počátku

Korelace lze v případě fitování polynomem výrazně omezit jednoduchým trikem (pokud je nezávislá proměnná $x$ vzorkována rovnoměrně): posuneme-li počátek doprostřed intervalu hodnot $x$, všechny korelace koeficientů u sudé a liché mocniny budou blízké nule. Je totiž $\sum_i (x_i-x_m)^{2k+1}=0$ při sčítání přes interval symetrický kolem $x_m$. Za těchto podmínek jsou liché a sudé sloupce modelové matice navzájem otrogonální.

In [19]:
xb=x-x.mean()
amatb=np.array([np.ones_like(xb),xb,xb**2]).T
hessb=amatb.T@amatb
par2b=np.linalg.inv(hessb)@amatb.T@j
par2b
Out[19]:
array([ 4.34487150e+02,  1.59308362e+01, -2.38671283e-01])
In [11]:
covb=np.linalg.inv(hessb)
errb=np.sqrt(covb.diagonal())
corb=covb/errb[np.newaxis,:]/errb[:,np.newaxis]
corb
Out[11]:
array([[ 1.        ,  0.        , -0.74565194],
       [ 0.        ,  1.        ,  0.        ],
       [-0.74565194, -0.        ,  1.        ]])

Omezení korelací způsobí zkrácení elips "ekvipotenciálních" kontur minimalizované funkce (viz demonstrace u fitování Drudeho modelu), což výrazně zmenší i velikost nejistot (hodnoty na diagonále kovarianční matice):

In [17]:
print("nejistoty bez posunuti: ",err)
print("nejistoty po posunuti:  ",errb)

nejistoty bez posunuti:  [4.32843141 0.74592287 0.0211732 ]
nejistoty po posunuti:   [0.23437702 0.01319909 0.00124831]

Obě modelové křivky (a tedy i hodnoty reziduí či jejich sumy) jsou samozřejmě identické - v obou případech fitujeme parabolu za stejných podmínek (stejné minimalizované veličiny). Rozdíl je na úrovni zaokrouhlovací chyby:

In [21]:
model2=amat@par2
model2b=amatb@par2b
max(abs(model2-model2b))
Out[21]:
1.7053025658242404e-12
In [29]:
xb=x-x.mean()
lig_ml=lambda p:-(j*np.log(abs(np.polyval(p[::-1],xb)))).sum()+np.polyval(p[::-1],xb).sum()
par1b=op.fmin(lig_ml,par2b)
par1b

Optimization terminated successfully.
         Current function value: -84309.680558
         Iterations: 84
         Function evaluations: 154
Out[29]:
array([ 4.26104370e+02,  1.58652646e+01, -1.78794424e-01])
In [23]:
hess3b=amatb.T@wmat@amatb
par3b=np.linalg.inv(hess3b)@amatb.T@wmat@j
par3b
Out[23]:
array([ 4.26616891e+02,  1.59395427e+01, -1.88620135e-01])

Pokud započítáme "poissonovské" váhy, korelace jsou o poznání větší - ortogonalita už není splněna při aplikaci wmat. Jednotlivé sloupce modelové matice by se daly ortogonalizovat obvyklými matematickými postupy s případným započtením váhy (je známo hodně ortogonálních sad polynomů), ve výsledku ale bude průběh modelové funkce stejný (a určované parametry ztratí svůj zjevný smysl).

In [34]:
cov3b=np.linalg.inv(hess3b)
err3b=np.sqrt(cov3b.diagonal())
cor3b=cov3b/err3b[np.newaxis,:]/err3b[:,np.newaxis]
cor3b
Out[34]:
array([[ 1.        , -0.00562375, -0.71434362],
       [-0.00562375,  1.        ,  0.52321168],
       [-0.71434362,  0.52321168,  1.        ]])