Recommend using Xarray/InferenceData

[hyunjimoon]

• "random variable as datatype" (Andrew's suggestion being implemented in Stan as rvars)

prior_pred_obs.mean(dim = "chain")

<xarray.Dataset>
Dimensions:                            (draw: 1000,
                                        production_start_rate_stocked_obs_dim_0: 20,
                                        production_rate_stocked_obs_dim_0: 20)
Coordinates:
  * draw                               (draw) int64 0 1 2 3 ... 996 997 998 999
Dimensions without coordinates: production_start_rate_stocked_obs_dim_0,
                                production_rate_stocked_obs_dim_0
Data variables:
    production_start_rate_stocked_obs  (draw, production_start_rate_stocked_obs_dim_0) float64 ...
    production_rate_stocked_obs        (draw, production_rate_stocked_obs_dim_0) float64 ...

prior_pred_obs.mean(dim = "draw")

<xarray.Dataset>
Dimensions:                            (chain: 1,
                                        production_start_rate_stocked_obs_dim_0: 20,
                                        production_rate_stocked_obs_dim_0: 20)
Coordinates:
  * chain                              (chain) int64 1
Dimensions without coordinates: production_start_rate_stocked_obs_dim_0,
                                production_rate_stocked_obs_dim_0
Data variables:
    production_start_rate_stocked_obs  (chain, production_start_rate_stocked_obs_dim_0) float64 ...
    production_rate_stocked_obs        (chain, production_rate_stocked_obs_dim_0) float64 ...

[hyunjimoon]

prior_pred_obs['production_start_rate_stocked_obs'].mean('draw')

<xarray.DataArray 'production_start_rate_stocked_obs' (chain: 1,
                                                       production_start_rate_stocked_obs_dim_0: 20)>
array([[-4.40146950e-03,  3.36716915e-03,  7.66082443e-03,
         1.60187649e-02, -3.77674004e-03,  4.27628952e-03,
         1.13146864e-03, -4.19492450e-03, -6.59530433e-03,
        -4.86387919e-03,  5.67682390e-03,  5.63749278e+09,
         5.10863331e+09,  2.47952251e+09,  4.95773346e+06,
         5.53670186e-01, -1.44968518e-03,  1.15607441e-02,
         5.49195196e-03,  7.08743052e-03]])
Coordinates:
  * chain    (chain) int64 1
Dimensions without coordinates: production_start_rate_stocked_obs_dim_0

[hyunjimoon]

Questions to be discussed with @Dashadower (in blue or ?)

Summary

why

• don't want to reset_index pandas (Q1.)
• labeled array (Q2. relevant to stan's ragged array)

what

multi-dim labeled array which would allow temporal and cross-sectional comparison

how

with two components (DataArray, DataSet) each with four (partially) hierarchical servants

• values: a numpy.ndarray holding the array’s values
• dims: dimension names for each axis e.g., ('x', 'y', 'z') (chain, draws)
• coords: dict-like container of arrays (coordinates) that label each point e.g., 1-dimensional arrays of numbers, datetime objects or strings
• attribute: OrderedDict to hold arbitrary metadata (attributes), metadata
  Q3. why does dims use set (), coords use dict-like container, attribute use OrderedDict?

A dataset is simply an object containing multiple DataArrays indexed by variable name

Let's start with some very basic examples to see it in action. We then proceed to a more realsistic example.
e.g.

coordinates and data variables; the latter is equivalent to pandas' column name. arviz (Bayes computation visualization library) used netCDF data structure as it need labels (posterior, prior, prior predictive, loglik) for each of the parameters each of which has 1,000 samples. coordinate,

From xarray cheat sheet

Q4. what is the difference btw point coordinate and single coordinate dimension

Q5. So, don't use loc when we are not changing the original value is their message?
The.loc[] method is another way to slice data, but I advise against using it unless you're actually changing values for a particular slice of the data, as it does not return a deep copy like the.sel method does, so it still points to the original data (unless you save the slice to another variable, in which case, still don't risk it!!!).

A shallow copy constructs a new compound object and then (to the extent possible) inserts references into it to the objects found in the original.
A deep copy constructs a new compound object and then, recursively, inserts copies into it of the objects found in the original.

Q5. expand dimension VS broadcasting

Q6.
A typical use case is when you want to keep several related data variables in one data structure but when they have inconsistent dimensions. For example, you want to store SDI and mortality in a dataset together, but mortality has age-group and sex dimensions, while SDI does not. They do however share year and location as dimensions.

Merging two or more dataarrays into one dataset: Note that dataarrays have to be named

[hyunjimoon]

In issue #17, can we store both param_classify and data_draws2data in one dataArray or dataSet?

# da = xr.DataArray(
#     data=temperature,
#     dims=["time", "classification"],
#     coords=dict(
#         lon=(["x", "y"], lon),
#         lat=(["x", "y"], lat),
#         time=time,
#         classification = param_classify,
#         reference_time=reference_time,
#     ),
#     attrs=dict(
#         description="mngInven.",
#     ),
# )

param_classify = {
    "est_params" : ("inventory_adjustment_time", "minimum_order_processing_time"),
    "ass_param_lst" : ("customer_order_rate", "inventory_coverage", "manufacturing_cycle_time", "safety_stock_coverage", "time_to_average_order_rate", "wip_adjustment_time"),
    "obs_stock_lst" : ("production_start_rate_stocked", "production_rate_stocked"),
    "driving_data_lst" : ("customer_order_rate", "ran_norm1", "ran_norm2"),
}

data_draws2data = {
    "n_t": n_t,
    "production_start_rate_stocked_obs": placeholder,
    "production_rate_stocked_obs":  placeholder,
    "customer_order_rate": customer_order_rate.loc[:n_t, 'y'].values,
    "standard_normal_for_process_noise": ran_norm1,
    "beta22shift1_4_production_rate_m_noise": ran_beta1,
    "beta22shift1_4_production_start_rate_m_noise": ran_beta2,
}

[hyunjimoon]

The shape of initial value is # of prior draws (S = 30) whereas the shape of stock vector is # of prior draws (S = 30)* # of time_window (n_t = N). Coordinates are chain, draw and each variable names are treaded as dims without coords.

<xarray.DataArray 'production_rate_stocked__init' (chain: 1, draw: 10)>
array([[198624., 198624., 198624., 198624., 198624., 198624., 198624.,
        198624., 198624., 198624.]])
Coordinates:
  * chain    (chain) int64 1
  * draw     (draw) int64 0 1 2 3 4 5 6 7 8 9
ds.production_rate_stocked
Out[26]: 
<xarray.DataArray 'production_rate_stocked' (chain: 1, draw: 10,
                                             production_rate_stocked_dim_0: 30)>
array([[[1.11186e+04, 9.58473e+03, ..., 7.70323e+08, 5.39617e+08],
        [1.11180e+04, 9.58418e+03, ..., 7.60463e+08, 5.32709e+08],
        ...,
        [1.11196e+04, 9.58553e+03, ..., 7.77838e+08, 5.44881e+08],
        [1.11194e+04, 9.58534e+03, ..., 7.58094e+08, 5.31050e+08]]])
Coordinates:
  * chain    (chain) int64 1
  * draw     (draw) int64 0 1 2 3 4 5 6 7 8 9
Dimensions without coordinates: production_rate_stocked_dim_0

ds.dims

Frozen({'chain': 1, 'draw': 10, 'initial_outcome_dim_0': 7, 'integrated_result_dim_0': 30, 'integrated_result_dim_1': 7, 'inventory_dim_0': 30, 'production_rate_stocked_dim_0': 30, 'work_in_process_inventory_dim_0': 30, 'expected_order_rate_dim_0': 30, 'production_start_rate_stocked_dim_0': 30, 'backlog_dim_0': 30, 'process_noise_dim_0': 30, 'production_rate_stocked_obs_dim_0': 30, 'production_start_rate_stocked_obs_dim_0': 30})

ds.coords

  * chain    (chain) int64 1
  * draw     (draw) int64 0 1 2 3 4 5 6 7 8 9

Note 'initial_outcome_dim_0': 7, meaning ragged array; not every vectors are indexed along time axis and each dims has its own index. This explains why to_dataframe automatically creates multi-index in pandas i.e. cartesian product of every chosen dimension. For instance, when production_rate_stocked_obs_dim_0 and production_start_rate_stocked_obs_dim_0 are chosen, the returned dataframe's length (shape[0]) is 900.

[hyunjimoon]

Arviz's datatype explains this with: posterior, posterior_predictive, sample_stats, prior, observed_data.
https://arviz-devs.github.io/arviz/getting_started/XarrayforArviZ.html#an-introduction-to-each

chain and draws are default coordinates, but each parameter dimension is dimension without coordinates which can help expanding hierarchy.

[hyunjimoon]

This can be connected to visual plots provided by arviz and online interactive exploration tool shinystan.

[hyunjimoon]

@tomfid may I carefully ask whether vensim can support inferencedata as its output or not? Also, if thankfully so, around when? I am not asking for quick answer, but please let me introduce a feedback loop we can operationalize in November.

Vensim's commitment is centripetal force which will allow me to capture great talents around me. Honestly, I need great help in developing (e.g. stanify issues) and November where hierarchical Bayesian and simulation-based test will be covered in @hazhirr's seminar is once in a year opportunity to launch a platform. 24 people who committed their time for this topic can help testing a platform on which system dynamics and Bayesian modeling are connected with standardized inferencedata. This platform can realize actionable workflow which can benefit Vensim as well.

Modularized development plans "Iet's develop plots given inferencedata (which Vensim will give)" is my persuasion strategy which has 100% success rate so far! as @cdigenna and @RossFW thankfully offered help as they appreciated the value of visualization for model critique, evaluation, and comparison. We aim to start designing our code from assumed inferencedata datastructure using az.load_arviz_data command. However, I am being a fraud if the assumption (given inferencedata) is wrong, or making it right is too far future. As SD modelers are experts in feedback loop of data-informed decision-making so I can't do this alone. Not only plot, but also how we use the plots need deep discussion.

I've gathered some info that can be of reference to our team :)

General

SciPy 2019 conference talk is one of the best I've seen which introduces using arviz for Bayesian testings. It

• motivates probabilistic programming with golf example
• compares Bayes with ML
• trains us to classify MCMC traceplot in to big step, small step, good (which is HMC).
• helps vensim-stan connection with code-based details of infercedatatype and xarray
• mesmerize me with raining visualization which can go well with sequential updates!

Vensim-Stan connection team

@hyunjimoon @tomfid

• Arviz's NetCDF connection support
• Arviz schema

InferenceData Plot team

@hyunjimoon @cdigenna @RossFW

• we can start from Arviz's example which starts from pre-loaded "centered_eight" dataset which includes prior, observed data, posterior samples for hierarchical Bayesian model

15.879 seminar team

@hyunjimoon

• Encouraging 24 testers (e.g. brief inclusion in assignment 10), and I would need Hazhir's help on this.

[hyunjimoon]

In fact, it was pysd developer @enekomartinmartinez who recommended xarray (familty of netCDF and inferencedata). Pysd supports netCDF with this function.

@jandraor and I planned to connect inferencedata with SBC package where we can use coverage plots. Detailed plans are documented in this SBC package issue.

[hyunjimoon]

dims is axis name, coords each data var's non-null dims

[hyunjimoon]

Questions to discuss with @OriolAbril on 12W2

• when to use xr.concat([])
• how to plot by time?
• the difference between grouping dataarrays into dataset VS grouping dataset into groups of dataset
• reason why groups of dataarray (i.e. dataset) is not enough for inferencedata
• could we infer hier_est_param_name from .nc?
• In plot_qoi, how to replace to_dataframe?
• do you agree with bayesflow's implementation on N/B > 20 when B = L/2? I have four chains, with L = 250 (or 1000?)

[hyunjimoon]

• independence btw dim and coord different group (chain, draw should be independent) - between prior and posterior. posterior and posterior predictive (too big)
• same result isel, sel (coordinate is defunct)
• xr.apply_ufunc complicated numpy function (linear algebra objects) to xarray object without losing dimension info

[hyunjimoon]

DecW3 Oriol, Angie

goal resetting

• plans to connect inferencedata usecase with simulation community (dynamic modeling, decision-based optimization) e.g. vensim (compartment dynamic simulation) and entity (agent-based modeling)
• detailed plans for three datasets here
• how to design output format (e.g. kwargs)?
• could inferencedata host the following three dataset examples with additional group concept?

action items

1. fed.agency data (theme: hierarchical modeling, five-layer)

• how to do multi-groupby?

2. imdb data (theme: agent-based modeling: title, year, person, character)

• when groups are needed?
• is kwargs dims?

indexed by (title, year). Interestingly this reminds me of stacking chain and draw into prior_draw.

3. sbc data (theme: verification and plots for dynamic simulation)

• why does sbc.prior doesn't have prior_draw as index?
• sbc next step usingwrapper + connecting with rvars

basic xarray questions

• where index is set from the following example?

dims = ('example', 'x', 'y')
combined = xarray.Dataset(
	data_vars = {'var': (dims, np.zeros(len(examples), 100, 200))}
	coords = {'example': examples}
	)
for example in examples:
	combined.loc[dict(example = example)] = create_an_xarray_dataset(example)	

• dims without coords

"Unindexed dimensions" is simply dimensions that don't have an associated coordinate variable. A better section title might be "Dimensions without coordinates". There's nothing inherently wrong with this: sometimes there aren't meaningful coordinates values for a dimension, and xarray handles this just fine.

• could we have coordinate that aren't from dims?
• multiple groupby agency.stack(c_new= [c1, c2], create_index = False).groupby('c_new`).count()
• 'c_new' or c_new (is named label string of variable name)
• how index system works under the hood (does xarray use pandas/numpy idx?) + does xarray use np.where where?
• why can we use sel only with indexes? (not year which is a coordinate indexed by titleId)
• when do we use swap_dims, "one option is to make year an index, then use sel(year=2000) (you can do that with swap_dims"
• Osvaldo's comment on bokeh's advantage (propagating subplots?)
• what is rcParams?

[hyunjimoon]

Template for PTSCV analysis
parameter[P] -> Time[N], Space [R] // Componet[C] -> Value[V]

# non-spatial subscript
variant_list = ["Delta", 'Omicron', "BA5"]
vax_status_list = ['Vx', 'NVx']

# spatial subscript
state_names = ['Alabama', 'Alaska', 'Arizona', 'Arkansas', 'California', 'Colorado', 'Connecticut', 'Delaware', 'Florida', 'Georgia', 'Hawaii', 'Idaho', 'Illinois', 'Indiana', 'Iowa', 'Kansas', 'Kentucky', 'Louisiana', 'Maine', 'Maryland', 'Massachusetts', 'Michigan', 'Minnesota', 'Mississippi', 'Missouri', 'Montana', 'Nebraska', 'Nevada', 'New Hampshire', 'New Jersey', 'New Mexico', 'New York', 'North Carolina', 'North Dakota', 'Ohio', 'Oklahoma', 'Oregon', 'Pennsylvania', 'Rhode Island', 'South Carolina', 'South Dakota', 'Tennessee', 'Texas', 'Utah', 'Vermont', 'Virginia', 'Washington', 'West Virginia', 'Wisconsin', 'Wyoming', 'District of Columbia', 'American Samoa', 'Guam', 'Northern Mariana Islands', 'Puerto Rico', 'United States Minor Outlying Islands', 'U.S. Virgin Islands']

param_lst =[]
for name in variant_list:
    param_lst.append(f'Variant Accuity Multiplier[{name}]')
for name in vax_status_list:
    param_lst.append(f'Immunity Loss Time[{name}]')
curr_param_list = []
for name in state_to_code_dict.keys():
    curr_param_list.append(f"Weight on Reported Probability of Infection[{name}]")
param_lst = param_lst + curr_param_list

dv = {
    "coords": {
        'value': {"dims": ("value"), "data": ["marginal", "average"]},
        'time': {"dims": ("time"), "data": [i+1 for i in range(1100)]},
        'params': {"dims": ("params"), "data": param_lst},
    },
    "dims": {"value", "time", "params"},
    "data_vars":{
        'VV': {"dims": "value", "data":[1,2]},
        'VV_1': {"dims": ("params", "value"), "data": np.zeros(shape = (62, 2))},
        'VV_t': {"dims": ("time", "value"), "data": np.zeros(shape = (1100, 2))},
    },
   "attrs": {"title": "vaccine value disaggregated by parameter, time, space, component"}
    }
vds = xr.Dataset.from_dict(dv)

[hyunjimoon]

1W2 Oriol, Angie

sbc-inferencedata

• dynamic model; below is time series interface we are designing

set_prior("stock[t, hier_param] ~ normal(mu[t, hier_param], noise_scale)")

• likelihood by timeseries components + how and when to use arviz wrapper?

arviz dev on rvar?

• arviz's plan on rvar connection? I have meeting with Paul, Martin, Shinyoung on SBC further dev next Thursday

misc. q on plots

• assign vs .where - do you operate with index or does pandas, xarray, arviz's handy command take cares of it?

• df[df['3AGY'].isin(['VA', 'AR'])] vs df[df.3AGY.isin(['VA', 'AR'])] (doesn't work)

• could we denote recent date to have higher alpha?

• do you recommend having date as int(left) or string(right)?
  

• general aesthetic skills of plots (hue, style from seaborn scatter plots ) from vicky.ipynb
  

[hyunjimoon]

example of vensim subscript

birth rate=
	0.08
	~	fraction/Year
	~		|

aging[age]=
	Population[age] / average time in cohort[age]
	~	person/Year
	~		|

initial population[age]=
	30, 50, 50, 60, 50
	~	person
	~		|

average time in cohort[age]=
	2, 10, 7, 40, 20
	~	Year
	~		|

births=
	Population[middle] * birth rate
	~	person/Year
	~		|

Population[infant]= INTEG (
	births - aging[infant],
		initial population[infant]) ~~|
Population[all but youngest]= INTEG (
	aging[previous cohort] - aging[all but youngest],
		initial population[all but youngest])
	~	person
	~		|

age:
	infant, child, teen, middle, old
	~	index
	~		|

all but youngest:
	child, teen, middle, old
	~	index
	~	age cohorts except infant
	|

previous cohort:
	infant, child, teen, middle -> all but youngest
	~	index
	~		|

********************************************************
	.Control
********************************************************~
		Simulation Control Paramaters
	|

FINAL TIME  = 100
	~	Year
	~	The final time for the simulation.
	|

INITIAL TIME  = 0
	~	Year
	~	The initial time for the simulation.
	|

SAVEPER  = 1
	~	Year
	~	The frequency with which output is stored.
	|

TIME STEP  = 1
	~	Year
	~	The time step for the simulation.
	|

[hyunjimoon]

for time and space template

import xarray as xr
import numpy as np
# non-spatial subscript
variant_list = ["Delta", 'Omicron', "BA5"]
vax_status_list = ['Vx', 'NVx']

# spatial subscript
state_names = ['Alabama', 'Alaska', 'Arizona', 'Arkansas', 'California', 'Colorado', 'Connecticut', 'Delaware', 'Florida', 'Georgia', 'Hawaii', 'Idaho', 'Illinois', 'Indiana', 'Iowa', 'Kansas', 'Kentucky', 'Louisiana', 'Maine', 'Maryland', 'Massachusetts', 'Michigan', 'Minnesota', 'Mississippi', 'Missouri', 'Montana', 'Nebraska', 'Nevada', 'New Hampshire', 'New Jersey', 'New Mexico', 'New York', 'North Carolina', 'North Dakota', 'Ohio', 'Oklahoma', 'Oregon', 'Pennsylvania', 'Rhode Island', 'South Carolina', 'South Dakota', 'Tennessee', 'Texas', 'Utah', 'Vermont', 'Virginia', 'Washington', 'West Virginia', 'Wisconsin', 'Wyoming', 'District of Columbia', 'American Samoa', 'Guam', 'Northern Mariana Islands', 'Puerto Rico', 'United States Minor Outlying Islands', 'U.S. Virgin Islands']

param_lst =[] #['Sensitivity to Weather']
curr_param_list = []
for name in state_names:
    curr_param_list.append(f"Weight on Reported Probability of Infection[{name}]")
param_lst.append(curr_param_list)
for name in variant_list:
    param_lst.append([f'Variant Accuity Multiplier[{name}]'])
for name in vax_status_list:
    param_lst.append([f'Immunity Loss Time[{name}]'])

dv = {
    "coords": {
        'value': {"dims": ("value"), "data": ["marginal", "average"]},
        'component': {"dims": ("component"),  "data": ['death', 'gdp', 'hospitalization']},
        'time': {"dims": ("time"), "data": [i+1 for i in range(1100)]},
        'params': {"dims": ("params"), "data": param_lst},
        'space': {"dims" : ("space"), "data": state_names}
    },
    "dims": {"value", "component", "time", "params"},
    "data_vars":{ #5+ 3
        'VV': {"dims": "value", "data":[1,2]},
        'VV_cv': {"dims": ("component", "value"), "data": np.zeros(shape = (3, 2))},
        'VV_sv': {"dims": ("space", "value"), "data": np.zeros(shape = (len(state_names), 2))},
        'VV_tv': {"dims": ("time", "value"), "data": np.zeros(shape = (1100, 2))},
        'VV_ctv': {"dims": ("component", "time", "value"), "data": np.zeros(shape = (3, 1100, 2))},
        'VV_1_pv': {"dims": ("params", "value"), "data": np.zeros(shape = (len(param_lst), 2))},
        'VV_1_pcv': {"dims": ("params", "component", "value"), "data": np.zeros(shape = (len(param_lst),3, 2))},
        'VV_1_ptv': {"dims": ("params", "time", "value"), "data": np.zeros(shape = (len(param_lst),1100, 2))},
        'epsilon_pv': {"dims": ("params", "value"), "data": np.zeros(shape = (len(param_lst), 2))},
        'epsilon_pcv': {"dims": ("params", "component", "value"), "data": np.zeros(shape = (len(param_lst), 3, 2))},
    },
   "attrs": {"title": "vaccine value disaggregated by parameter, time, space, component"}
    }
vds = xr.Dataset.from_dict(dv)

[hyunjimoon]

07.11 Oriol, Angie
1d arrays as things
coord as two dims: xarray (geo)

radon example (100 x 1000: radon concentration; chain, )
three 1d coordinate county, state, floor

unique observation: scalar

data_vars (actual parametr, posterior predictive samples) vs coordinates are indexing

data_vars > coords (linear reg: two parameters (alpha, beta, sigma), and one ; coord (chain, draws, obs_id))

posterior:
observation: obs_id

rugby - theme, county coordinate for posterior; observation would also have county coordinate (even if dgp is different by)-85 in 10000 obs (different);

tidy data (any data into df, specific obsrevation, column for obs- the rest is dim, coord, )
1d var (a,b,sig; four chains, thoundand draws, ) 41000 2d array (4000 single; chain, draw, a, b, ); 4100080 (3d)
41000*80 vs 320k tidyverse-constrainted (plot) assume data has tidy format (useful for plot, ragged from different sources; posterior, parmaeter in covariance (determinat, multiplied; flat data; dim1 and 2 into matrix))

data don't match tidy constraints (number of rows explode quickly; label coordinate, muliple index, ten times memory;), xarray is less plot focused at first

reducing the space dimension ()

data - how can i generate; add new coordinate for plotting (reshaping - in xarray can't reduce; plot posterior means (reduce chain and draw ))

pre-emptively create variables data_vars (numpy - initialize; feel that array and pre-emptive create array); the morelooping
looping
after initialize (dataset square bracket[variable] = data)

.loc :slice of the variable

[hyunjimoon]

hyunjimoon/VaccineMisinf#15