Covariance analysis and estimators
[1]:
import lenapy
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
/work/scratch/env/fourests/.conda/envs/lenapy/lib/python3.11/site-packages/esmpy/interface/loadESMF.py:92: VersionWarning: ESMF installation version 8.6.1, ESMPy version 8.6.0
warnings.warn("ESMF installation version {}, ESMPy version {}".format(
Open the global mean sea level timesseries
[2]:
gmsl=xr.open_dataset('../../data/MSL_wo_seasonal_signal.nc',engine='lenapyNetcdf')
data=gmsl.msl
Build the covariance matrix corresponding to a white noise between to dates
[3]:
a=data.lntime.covariance_analysis()
a.add_errors('random',1,None,lenapy.utils.time.JJ_to_date(19000),lenapy.utils.time.JJ_to_date(24000),None)
a.visu()
Build the covariance matrix corresponding to bias uncertainty at a given date
[4]:
a=data.lntime.covariance_analysis()
a.add_errors('bias',1,lenapy.utils.time.JJ_to_date(19116),None,None,None,'centered')
a.visu()
Build the covariance matrix corresponding to a drift uncertainty between to dates
[5]:
a=data.lntime.covariance_analysis()
a.add_errors('drift',0.0001,None,lenapy.utils.time.JJ_to_date(19000),lenapy.utils.time.JJ_to_date(24000),None,'right')
a.visu()
Build the covariance matrix corresponding to a correlated noise uncertainty over 2 years (73 days) between to dates
[6]:
a=data.lntime.covariance_analysis()
a.add_errors('noise',1,None,lenapy.utils.time.JJ_to_date(19000),lenapy.utils.time.JJ_to_date(24000),np.timedelta64(730,'D'),None)
a.visu()
Build a covariance matrix from a yaml file
[7]:
error_prescription="../../data/errors.yaml"
cov=data.lntime.covariance_analysis()
cov.read_yaml(error_prescription)
cov.visu()
Perform OLS regressions on the gmsl timeseries
[12]:
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,5))
data.plot(ax=ax1,label='GMSL')
# Degree 1
est=data.lntime.OLS(degree=1)
est.estimate.plot(ax=ax1,label='degree=1')
est.residuals.plot(ax=ax2,label='degree=1')
# Degree 2
est=data.lntime.OLS(degree=2)
est.estimate.plot(ax=ax1,label='degree=2')
est.residuals.plot(ax=ax2,label='degree=2')
ax1.legend()
ax2.legend()
print(est.uncertainties)
<xarray.DataArray (degree: 3)> Size: 24B
array([4.44151369e-02, 1.04455671e-10, 4.12706472e-19])
Coordinates:
* degree (degree) int64 24B 0 1 2
Perform an GLS regressions on the gmsl timeseries with the built covariance matrix
[13]:
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,5))
data.plot(ax=ax1,label='GMSL')
# Degree 1
est=data.lntime.GLS(degree=1,sigma=cov.sigma)
est.estimate.plot(ax=ax1,label='degree=1')
est.residuals.plot(ax=ax2,label='degree=1')
# Degree 2
est=data.lntime.GLS(degree=2,sigma=cov.sigma)
est.estimate.plot(ax=ax1,label='degree=2')
est.residuals.plot(ax=ax2,label='degree=2')
ax1.legend()
ax2.legend()
print(est.uncertainties)
<xarray.DataArray (degree: 3)> Size: 24B
array([9.33275009e-04, 5.46918674e-12, 1.16288723e-20])
Coordinates:
* degree (degree) int64 24B 0 1 2
[ ]: