<jupyter_start><jupyter_text>MIC Demo 1 - Basic steps for measurementThis simple demonstration of the MIC toolbox uses two simulated bivariate VAR(2) models from the ["Macroeconomic simulation comparison with a multivariate extension of the Markov Information Criterion"](https://www.kent.ac.uk/economics/documents/research/papers/2019/1908.pdf) paper. These are the first two settings from the VAR validation exercises. The two simulated datasets are located in `data/model_1.txt` and `data/model_2.txt`. In addition this, one of these two models has been used to generated an 'empirical' dataset `data/emp_dat.txt`. The purpose of the demonstration is to show how to run the MIC and see if we can figure out which of models 1 or 2 is the true model.The purpose of this first part is to outline the individual steps required to obtain a MIC measurement on a single variable in a multivariate system. As a full multivariate measurment requires several runs of the algorithms, this is best done in parallel, which will be covered in the second demonstration notebook.We start with the setup, including the toolbox import:<jupyter_code>import time
import numpy as np
from scipy.stats import pearsonr
import mic.toolbox as mt<jupyter_output><empty_output><jupyter_text>Stage 0 - Discretising the dataThe first task that needs to be done is to discretise the two variables in the system (denoted $x^1$ and $x^2$). In order to do so, we need to provide the following information:- `lb` and `ub` : Bounds to the range of variation. - `r_vec` : Binary discretisation resolution of the variables- `hp_bit_vec` : High priority bits - Number of bits to prioritise in the permutationWe can then call the binary quantisation function in the toolbox, `mt.bin_quant()`, and look at the result of the discretisation diagnostics to ensure that settings above are chosen so that the discretisation error is i.i.d uniformly distributed.<jupyter_code>lb = [-10,-10]
ub = [ 10, 10]
r_vec = [7,7]
hp_bit_vec = [3,3]

# Load 'empirical' data 
path = 'data/emp_data.txt'
emp_data = np.loadtxt(path, delimiter="\t") 

# Pick first replication (columns 1 and 2) - just as an example
dat = emp_data[:,0:2]

# Run the discretisation tests (displays are active by passing any string other than 'notests' or 'nodisplay')
data_struct_emp = mt.bin_quant(dat,lb,ub,r_vec,'')

# Check the correlation of the high-priority bits (example of 'notests' here)
data_struct_hp = mt.bin_quant(dat,lb,ub,hp_bit_vec,'notests')
dat_bin = data_struct_hp['binary_data']
hp_dat = dat - np.asarray(dat_bin.errors)

print('\n Correlation of high priority bits with raw data:')
for n in range(2):
    corr = pearsonr(dat[:,n],hp_dat[:,n])
    print(' Variable {:1d}: {:7.4f}'.format(n+1,corr[0]))<jupyter_output>+-------------------------------------------------+
|            Quantization diagnostics             |
+-------------------------------------------------+
 N° of observations:       1000

 Var N°:             1       2
 Lower bound:    -10.000 -10.000
 Upper bound:     10.000  10.000
 Min obs.:        -5.259  -4.709
 Max obs.:         6.160   4.381
 Out of bounds:        0       0
 Resolution:           7       7
 Standard dev.:    1.737   1.580
 Discr. unit:      0.156   0.156

 Signal-to-Noise ratio for quantisation (in dB):
 Var N°:             1       2
 Theoretical:     31.710  30.886
 Effective:       25.748  24.540

+-------------------------------------------------+
 Kolmogorov-Smirnov test for uniformity
 H0: Distribution of errors is uniform

 Var N°:             1       2
 KS statistic:     0.016   0.059
 P-value:          0.999   0.059

+-------------------------------------------------+
 Ljung-Box tests for autocorrelation of errors
 H0: The errors are independently distr[...]<jupyter_text>We can see here that under the settings picked above, the quantisation errors for both series are indeed uniformly distributed (KS test not rejected), independent (LB test is not rejected) and the errors are not correlated with the discretisation levels. Furthermore, a discretisation using only the first three bits of each variable (the 'high priority bits') already has a 98% correlation with the raw data. This suggests that the settings are appropriate for the analysis. Stage 1 - learning model probabilities from the dataThe important parameters to choose at this stage relate to the size of the tree that we want to build with the simulated data. The parameters of interest are:- `mem` : the maximum amount of nodes we can initialise. As trees tend to be sparse, this does not have to match the theoretical number of the trees $2^D$.- `d`: maximum depth of the context trees (in bits). Combined with `mem`, these implement a cap on the amount of the amount of memory that can be used.- `lags`: The number of lags in the Markov process being used to learn the probabilities.<jupyter_code>mem = 200000
d = 24
lags = 2<jupyter_output><empty_output><jupyter_text>Choosing 2 Markov lags and 14 bits per observations means that the full context will be 28 bits long (not counting contemporaneous correlations). Given a maximum tree depth $D=24$, it is clear that some of the context will be truncated. We therefore need to permute the bits in the context to prioritise the most important ones and ensure only the least informative ones get truncated. In order to do so, the next step is to generate a permutation the context bits. As stated above, we are only demonstrating a single run of the MIC algorithm, and we choose to predict the first variable conditional on the context and the current value of the second variable. This is declared via the `var_vec` list below. To clarify the syntax, the first entry in the `var_vec` list identifies the variable to predict (1 in this case), and any subsequent entries in the list indentify contemporaneous conditioning variables (2 in our case)This will allow us to determine the value of $\lambda^1 (x_t^1 \mid x_t^2, \Omega_t)$. It is important to note that to get the full MIC measurement, will need to run the algorithm again. The steps are essentially the same, and this will be covered in the second part of the demonstration.<jupyter_code>num_runs = 2
var_vec = [1,2]  # This is the critical input, it governs the conditioning order.

perm = mt.corr_perm(dat, r_vec, hp_bit_vec, var_vec, lags, d)<jupyter_output>+------------------------------------------------------+
 Predicting variable:   1
 Contemporaneous conditioning vars: 
	[2]
 Number of conditioning lags:   2


 --------- Iteration   1 ---------
 Select lag   1 of var   1
 Correlation:  0.7443
 --------- Iteration   2 ---------
 Select lag   2 of var   2
 Correlation:  0.4978
 --------- Iteration   3 ---------
 Select lag   2 of var   1
 Correlation:  0.4974
 --------- Iteration   4 ---------
 Select lag   1 of var   2
 Correlation:  0.3625
 --------- Iteration   5 ---------
 Select contemporaneous var   2
 Correlation:  0.0899
+------------------------------------------------------+<jupyter_text>We now have all the elements required to train the tree. For the purpose of the demonstration the two simulated data files contain two training sets of 125,000 observations for each variable $x^1$ and $x^2$. The first set is located in the first two columns of the traing data, while the second set is located in columns 3 and 4. This division into two training sets is done in order to illustrate:- How to initialise a new tree on the first series- How to update an existing tree with further dataAs stated above, we are only carrying out a single run here, so we choose to learn the probabilities for the 1st model only. Once again, to get a measurement for the second model will require further runs which we do in the second part of the demonstration.<jupyter_code># Load model data
path = 'data/model_1.txt'
sim_data = np.loadtxt(path, delimiter="\t") 

# Pick a tag for the tree (useful for indentifying the tree later on)
tag = 'Model 1'

# Discretise the training data. 
sim_dat1 = sim_data[:,0:2]
data_struct = mt.bin_quant(sim_dat1,lb,ub,r_vec,'notests') # Note the 'notests' option
data_bin = data_struct['binary_data']

# Initialise a tree and train it, trying to predict the 1st variable
var = var_vec[0]
output = mt.train(None, data_bin, mem, lags, d, var, tag, perm)<jupyter_output>+------------------------------------------------------+
|       Context Tree Weighting on training series      |
+------------------------------------------------------+
           Empty tree structure initialised
+------------------------------------------------------+
	 10% complete	   10.1220 seconds
	 20% complete	   18.8578 seconds
	 30% complete	   27.1562 seconds
	 40% complete	   35.4127 seconds
	 50% complete	   43.4386 seconds
	 60% complete	   51.9838 seconds
	 70% complete	   59.8614 seconds
	 80% complete	   68.0705 seconds
	 90% complete	   75.2379 seconds
	100% complete	   82.3344 seconds
+------------------------------------------------------+
 Total hashing time required      :    59.1277
 Average per context              : 4.7303e-04
 Total CTW updating time required :    19.7348
 Average per context              : 1.5788e-04
+------------------------------------------------------+<jupyter_text>Let's now update the tree with the second run of training data to see the difference in syntax and output<jupyter_code># Discretise the second run of training data
sim_dat1 = sim_data[:,2:4]
data_struct = mt.bin_quant(sim_dat1,lb,ub,r_vec,'notests') # Note, we are not running discretisation tests
data_bin = data_struct['binary_data']

# Extract the tree from the previous output and train it again. Only the 1st argument changes
T = output['T']
output = mt.train(T, data_bin, mem, lags, d, var, tag, perm)<jupyter_output>+------------------------------------------------------+
|       Context Tree Weighting on training series      |
+------------------------------------------------------+
               Using "Model 1" structure
+------------------------------------------------------+
	 10% complete	    7.5647 seconds
	 20% complete	   15.5909 seconds
	 30% complete	   23.4861 seconds
	 40% complete	   31.5217 seconds
	 50% complete	   39.5597 seconds
	 60% complete	   47.5197 seconds
	 70% complete	   55.6947 seconds
	 80% complete	   63.7563 seconds
	 90% complete	   71.4348 seconds
	100% complete	   78.3093 seconds
+------------------------------------------------------+
 Total hashing time required      :    54.0317
 Average per context              : 4.3226e-04
 Total CTW updating time required :    20.4773
 Average per context              : 1.6382e-04
+------------------------------------------------------+<jupyter_text>Notice how the header of the output has changed, using the tag to flag that we are updating an existing tree. We are done training the tree, let's get some descriptive statistics.<jupyter_code># Use the built-in descriptive statistic method to get some diagnostics
T = output['T']
T.desc()<jupyter_output>+------------------------------------------------------+
               Tree diagnostics: Model 1
+------------------------------------------------------+
 Variable 1 of 2
 Training observations:           249996
+------------------------------------------------------+
 Memory locations allocated:      200000 nodes
 Memory usage per node:             2040 Bytes
 Memory space required:        408000000 Bytes
 Effective memory useage:      147628680 Bytes
+------------------------------------------------------+
 Resolution:             7 bits
 Tree depth:            24 bits
 Leaf nodes:         16724 nodes
 Branch nodes:       55643 nodes
 Free nodes:        127633 nodes
+------------------------------------------------------+
 N° of failed allocations:           0
 Maximum N° of leaf tries:           7
 Average N° of leaf tries:      1.3716
 Maximum N° of branch tries:        10
 Average N° of branch tries:    1.3506
+------------------------------------------------------+
 Number of [...]<jupyter_text>We can see that the tree has used about 1/3 of the initial node allocation, so we have plenty of margin on memory. This will typically change if more variables are included. There is an element of trial and error to figuring out how much memory to allocated, which is why this diagnostic is useful. It is important to note that the algorithm can cope with failed node allocations (situations where the algorithm attempts to allocate a node to memory but fails), as it has a heuristic that allows it to 'skip' failed nodes, at the cost of introducing an error in the probabilities. Furthermore, because the tree implements a pruning and rollout mechanism, nodes are only allocated when they are not on a single branch path. This means that node allocation failures due to lack of memory will typically occur for very rare eventsBoth of these mechanisms mean that memory allocation failure is graceful and the odd failure will not impair the measurement. It is nevertheless a good idea to check `T.desc()` to ensure that failed allocations are not too numerous. Stage 2 - Scoring the empirical data with the model probabilitiesWe have already loaded the empirical data when running the discretisation tests. For the purposes of this demonstration, we have 10 replications of 1000 observations each, in order to show the consistency of the measurement. We will therefore loop the steps required to score these series over the 10 replications. In a normal application with only one emprical dataset, this loop is of course not needed!- Discretise the empirical data into `data_stuct_emp`- Extract the binary data from the dictionary - Pass the binary data alongside the tree `T` to the score function. The function knows which variable to score as this is given in the tree.- Correct the score by removing the estimate of the bias (measured using the Rissanen bound correction).<jupyter_code>scores = np.zeros([998,10])

for j in range(10):
    loop_t = time.time()

    # Discretise the data
    k = 2*j
    dat = emp_data[:,k:k+2]
    data_struct_emp = mt.bin_quant(dat,lb,ub,T.r_vec,'notests')
    data_bin_emp = data_struct_emp['binary_data']

    # Score the data using the tree
    score_struct = mt.score(T, data_bin_emp)

    # Correct the measurement
    scores[:,j] = score_struct['score'] - score_struct['bound_corr']

    print('Replication {:2d}: {:10.4f} secs.'.format(j,time.time() - loop_t))

flt_str = '    {:7.2f}'*10    
print('\n Scores obtained ')
print(flt_str.format(*np.sum(scores,0)))<jupyter_output>Replication  0:     0.7932 secs.
Replication  1:     0.8078 secs.
Replication  2:     0.8371 secs.
Replication  3:     0.8248 secs.
Replication  4:     0.8004 secs.
Replication  5:     0.8059 secs.
Replication  6:     0.8500 secs.
Replication  7:     0.8098 secs.
Replication  8:     0.8036 secs.
Replication  9:     0.7955 secs.

 Scores obtained 
    4793.87    4827.00    4794.54    4795.12    4784.41    4851.00    4788.13    4843.35    4805.88    4753.12
{'fraction_non_alphanumeric': 0.14562328904319957, 'fraction_numerical': 0.04473526073312412, 'mean_word_length': 2.4978981784212984, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 27, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '5060', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4373, 'n_tokens_neox': 3857, 'n_words': 1860, 'path': 'demo/MIC Demo 1 - Basic steps for measurement.ipynb', 'repo_name': 'Sylvain-Barde/mic-toolbox', 'size': '14977'}