Quickstart: the wrapper for the parameters estimator in Matlab

Expected time to run through: 5 mins

In this tutorial, we demonstrate how to run the Matlab code of the Bayesian estimator (see Tolwinski-Ward et al., 2013) for VS-Lite with PyVSL.

%load_ext autoreload
%autoreload 2

import PyVSL
import pandas as pd
import numpy as np

Load test data

First of all, we load the test data (can be downloaded with this link), which contains T and P at a grid point of the CCSM4 past1000 simulation. Note that the unit for T is degree C and that for P is accumulated monthly precipitation (mm/month).

data_dict = pd.read_pickle('./data/test_T_P.pkl')
time = data_dict['time']
T = data_dict['T'] # unit: degC
P = data_dict['P'] # unit: mm/month

print(time)
print(np.min(T), np.max(T))
print(np.min(P), np.max(P))

[ 850.04109589  850.12328767  850.20273973 ... 1850.7890411  1850.8739726
 1850.95616438]
-31.395203 16.369354
1.3947426 362.0332

Run the wrapper

We first generate pseudo-TRW and then use it to test the Bayesian estimator of the parameters.

res = PyVSL.VSL_M(
        850, 1850, 52.3,                  # the starting year and ending year of the input T & P, along with the latitude
        T, P,
        T1=3, T2=15, M1=0.05, M2=0.2,  # parameters of the thresholds
)

pseudoTRW = res['trw']
print(pseudoTRW)

[ 0.18525126  1.71342774  2.0978884  ... -1.01238715  0.08412163
 -0.98738839]

Now we estimate the paramters.

%%time

res_dict = PyVSL.est_params(T, P, 52.3, pseudoTRW)
print(res_dict)

Working on chain 1 out of 3...
Working on chain 2 out of 3...
Working on chain 3 out of 3...
warning: implicit conversion from numeric to char
warning: called from
    estimate_vslite_params_v2_3 at line 456 column 9
    _pyeval at line 57 column 30

warning: implicit conversion from numeric to char
warning: called from
    estimate_vslite_params_v2_3 at line 457 column 9
    _pyeval at line 57 column 30

warning: implicit conversion from numeric to char
warning: called from
    estimate_vslite_params_v2_3 at line 458 column 9
    _pyeval at line 57 column 30

warning: implicit conversion from numeric to char
warning: called from
    estimate_vslite_params_v2_3 at line 459 column 9
    _pyeval at line 57 column 30

warning: implicit conversion from numeric to char
warning: called from
    estimate_vslite_params_v2_3 at line 461 column 13
    _pyeval at line 57 column 30

warning: implicit conversion from numeric to char
warning: called from
    estimate_vslite_params_v2_3 at line 466 column 9
    _pyeval at line 57 column 30

    Rhat for T1, T2, M1, M2, sigma2rw:
    4.4953   85.2531    0.9990    0.9991   29.1848
Gelman and Rubin metric suggests MCMC has not yet converged to within desired threshold;
Parameter estimation code should be re-run using a greater number of MCMC iterations.
(See 'nsamp' advanced input option.)
{'T1': 2.9999204401190873, 'T2': 15.000326318028137, 'M1': 0.03217521381599234, 'M2': 0.2300052886329777}
CPU times: user 41.4 ms, sys: 21.8 ms, total: 63.2 ms
Wall time: 1min 16s