There have been proposals to use xarray as a common language for pymc3, pystan, and pymc4. This library might be a good place to start that by implementing utility functions for translating pymc3's Multitrace and pystan's OrderedDict into xarray objects, and then having all plotting functions work with xarrays.

Hi community, I have two models that I would like to compare using LOO and WAIC.

Model 1:

with pm.Model() as Model_1:
    α_v = pm.HalfCauchy('α_v',2)
    α_a = pm.HalfCauchy('α_a', 2) 
    α_max_ampl = pm.HalfCauchy('α_max_ampl',2)  
    α_tau = pm.HalfCauchy('α_tau', 2)     
    α_ter = pm.HalfCauchy('α_ter', 2)  
    
    β = pm.Normal('β', 0.334, 0.1)

    γ_v = pm.Normal('γ_v', params_group[-1,0], .1)
    γ_a = pm.Normal('γ_a', params_group[-1,1], .1) 
    γ_ter = pm.Normal('γ_ter', params_group[-1, 2], .1)     
    γ_max_ampl = pm.Normal('γ_max_ampl', params_group[-1,3], .1)
    γ_tau = pm.Normal('γ_tau', params_group[-1,4], .1)     

    μ_v = -α_v*pm.math.exp(-β*Age_years_group) + γ_v
    μ_a = α_a*pm.math.exp(-β*Age_years_group) + γ_a      
    μ_max_ampl = α_max_ampl*pm.math.exp(-β*Age_years_group) + γ_max_ampl
    μ_tau = α_tau*pm.math.exp(-β*Age_years_group) + γ_tau 
    μ_ter = α_ter*pm.math.exp(-β*Age_years_group) + γ_ter     
    
    σ_v = pm.HalfNormal('σ_v', .1)
    σ_a = pm.HalfNormal('σ_a', .1)
    σ_max_ampl = pm.HalfNormal('σ_max_ampl', .1)
    σ_tau = pm.HalfNormal('σ_tau', .1)      
    σ_ter = pm.HalfNormal('σ_ter', .1)     
    
    y_v = pm.Normal('y_v', μ_v, σ_v, observed=params_group[:, 0])    
    y_a = pm.Normal('y_a', μ_a, σ_a, observed=params_group[:, 1])     
    y_ter = pm.Normal('y_ter', μ_ter, σ_ter, observed=params_group[:, 2]) 
    y_max_ampl = pm.Normal('y_max_ampl', μ_max_ampl, σ_max_ampl, observed=params_group[:, 3])     
    y_tau = pm.Normal('y_tau', μ_tau, σ_tau, observed=params_group[:, 4])     

Model_2:

with pm.Model() as Model_2:
    α_v = pm.HalfCauchy('α_v',2)
    α_a = pm.HalfCauchy('α_a', 2) 
    α_max_ampl = pm.HalfCauchy('α_max_ampl',2)  
    α_tau = pm.HalfCauchy('α_tau', 2)     
    α_ter = pm.HalfCauchy('α_ter', 2)  
    
    β_v = pm.Normal('β_v', 0.334, 0.1)#Kail (1991)
    β_a = pm.Normal('β_a', 0.334, 0.1)#Kail (1991)
    β_max_ampl = pm.Normal('β_max_ampl', 0.334, 0.1)#Kail (1991)
    β_tau = pm.Normal('β_tau', 0.334, 0.1)#Kail (1991)
    β_ter = pm.Normal('β_ter', 0.334, 0.1)#Kail (1991)
      
    #asymptote (get the mean value of adults)
    γ_v = pm.Normal('γ_v', params_group[-1,0], .1)
    γ_a = pm.Normal('γ_a', params_group[-1,1], .1) 
    γ_ter = pm.Normal('γ_ter', params_group[-1, 2], .1)     
    γ_max_ampl = pm.Normal('γ_max_ampl', params_group[-1,3], .1)
    γ_tau = pm.Normal('γ_tau', params_group[-1,4], .1)     

    μ_v = -α_v*pm.math.exp(-β_v*Age_years_group) + γ_v
    μ_a = α_a*pm.math.exp(-β_a*Age_years_group) + γ_a      
    μ_max_ampl = α_max_ampl*pm.math.exp(-β_max_ampl*Age_years_group) + γ_max_ampl
    μ_tau = α_tau*pm.math.exp(-β_tau*Age_years_group) + γ_tau 
    μ_ter = α_ter*pm.math.exp(-β_ter*Age_years_group) + γ_ter     
    
    σ_v = pm.HalfNormal('σ_v', .1)
    σ_a = pm.HalfNormal('σ_a', .1)
    σ_max_ampl = pm.HalfNormal('σ_max_ampl', .1)
    σ_tau = pm.HalfNormal('σ_tau', .1)      
    σ_ter = pm.HalfNormal('σ_ter', .1)     
    
    y_v = pm.Normal('y_v', μ_v, σ_v, observed=params_group[:, 0])    
    y_a = pm.Normal('y_a', μ_a, σ_a, observed=params_group[:, 1])     
    y_ter = pm.Normal('y_ter', μ_ter, σ_ter, observed=params_group[:, 2]) 
    y_max_ampl = pm.Normal('y_max_ampl', μ_max_ampl, σ_max_ampl, observed=params_group[:, 3])     
    y_tau = pm.Normal('y_tau', μ_tau, σ_tau, observed=params_group[:, 4])           

I thus tried to compute:

compare_dict = {Model_1: trace, Model_2: trace_2}
compare_LOO = az.compare(compare_dict) 

and got the following error:

Traceback (most recent call last):

  File "<ipython-input-35-26e3f30f8929>", line 1, in <module>
    compare_LOO = az.compare(compare_dict)

  File "C:\Users\mservant\Anaconda3\lib\site-packages\arviz\stats\stats.py", line 211, in compare
    ics = ics.append([ic_func(dataset, pointwise=True, scale=scale)])

  File "C:\Users\mservant\Anaconda3\lib\site-packages\arviz\stats\stats.py", line 646, in loo
    log_likelihood = _get_log_likelihood(inference_data, var_name=var_name)

  File "C:\Users\mservant\Anaconda3\lib\site-packages\arviz\stats\stats_utils.py", line 412, in get_log_likelihood
    "Found several log likelihood arrays {}, var_name cannot be None".format(var_names)

TypeError: Found several log likelihood arrays ['y_v', 'y_a', 'y_ter', 'y_max_ampl', 'y_tau'], var_name cannot be None

Could you please help me solve this issue ? I am confused because I could compute LOO and WAIC using a previous version of arviz.

Description

Checklist

• [x] Follows official PR format
• [ ] Includes a sample plot to visually illustrate the changes (only for plot-related functions)
• [ ] New features are properly documented (with an example if appropriate)?
• [ ] Includes new or updated tests to cover the new feature
• [ ] Code style correct (follows pylint and black guidelines)
• [ ] Changes are listed in changelog

Added improved tooltip parameters to bokeh backend files of plots To fix #952

Description

This PR and #1747 summarizes my work done in Google Summer of Code 2021. I added 2 plots, plot_lm and plot_ts. The first one is covered in this PR and the next one in #1747. The project workflow was like this:

1. Started with a high-level design: Discussed the user API function with the mentors and other community members. Discussion includes input variable names, accepted datatypes, accepted values and other input details. Also showed a sample output visualization. Once it was approved, I moved on to opening this PR. This Github Gist shows the design and some discussions related to it. The design process took an initial week.

2. Submitted a prototype: This PR was earlier a prototype. This step is basically implementing the design decision made in the previous step. In addition to the user API, backend functions are also added. As ArviZ uses 2 backends, artist functions for both of them were added. This step took another week.

3. Improved it according to reviews: In this step, mentors review your code to provide feedback and improvement tips. Learned the best code practices in this step. Ideally, improving never ends. It is important to maintain the added code after it is merged to keep it bug-free. I aim to provide support after GSoC. This step was a bit lengthy and complex and thus took 2 weeks.

4. Tested the functionality: Added unit tests using pytest. Aimed and achieved to cover all of the added functionality under tests. Similar to step 2, added tests for both of the backends. Also, solved some simple to complex bugs that arise while testing. This step took another week.

5. Added examples: Examples are added to the docstring as well as to the examples folder. Please check out the files changed tab to know more about this step. This step was quick, took only half a week.

6. Added documentation: If you want to know how to use plot_lm, checkout out this blog. However, if you want to go on a low level to know the working in detail, I would suggest taking a look over the docstring in ArviZ docs and probably follow the comments sequentially. Another week was consumed in this step.

Output :

Checklist

• [x] Follows official PR format
• [x] Includes a sample plot to visually illustrate the changes (only for plot-related functions)
• [x] New features are properly documented (with an example if appropriate)?
• [x] Includes new or updated tests to cover the new feature
• [x] Code style correct (follows pylint and black guidelines)
• [ ] Changes are listed in changelog

Description

Start adding Dask compatibility starting from diagnostics/stats ufuncs.

Checklist

• [x] Follows official PR format
• [x] New features are properly documented (with an example if appropriate)?
• [ ] Includes new or updated tests to cover the new feature
• [ ] Code style correct (follows pylint and black guidelines)
• [ ] Changes are listed in changelog

Description

fixes #1053

Checklist

• [x] Does the PR follow official PR format?
• [x] Is the code style correct (follows pylint and black guidelines)?
• [x] Is the change listed in changelog?

Description

fixes #1801

Checklist

• [x] Follows official PR format
• [ ] Includes a sample plot to visually illustrate the changes (only for plot-related functions)
• [x] New features are properly documented (with an example if appropriate)?
• [ ] Includes new or updated tests to cover the new feature
• [x] Code style correct (follows pylint and black guidelines)
• [ ] Changes are listed in changelog

Changes Done: Replaced az.reloo with {func}~arviz.reloo Edited some grammatical mistakes Replaced PyStanSamplingWrapper terms with NumPyroSamplingWrapper Replaced "Stan Code" with "NumPyro Code" Replaced pystan_wrapper to numpyro_wrapper

Description

When running all tests locally ImportErrors are thrown if optional/extra dependencies are not installed. By skipping these tests when not installed, users can still run the entire test suite if desired. The entire test suite should still be run via CI.

This behavior is largely achieved by changing

import xx

to

xx = pytest.importorskip(xx)

which will try to import module xx, and will skip all tests in file if the module is not installed.

Individual tests may also be skipped using the decorator:

@pytest.mark.skipif(importlib.util.find_spec("xx") is None, reason="xx tests only required for CI")

The last remaining issue is numba, which is used in 20+ tests, used in a variety of files. One option is to make numba an entry in requirements-dev.txt instead of requirements-optional.txt, or add the skip decorators to each of the 20+ tests.

This PR addresses Issue #1112 .

Checklist

• [X] Follows official PR format
• [X] New features are properly documented (with an example if appropriate)?
• [X] Code style correct (follows pylint and black guidelines)
• [X] Changes are listed in changelog

~For now should solve only half of #220, the part related to constant_data.~ I would also like to support multiple log_likelihoods somehow (related to #771).

EDIT: constant data support was added in #814

EDIT

I update the description of what this PR does and what will still be pending.

Does

• Supports storage of multiple observed variables each with its own log likelihood data
  • dict
  • emcee
  • pymc
  • pystan
• Updates loo, waic... to get the log likelihood data from the new group. To do so, a function get_log_likelihood has been added in stats_utils. From here on, log likelihood values should be accessed using this helper function (it could even become public or a method of inference data, we have to think about this)

Note: lp variable is left in sample_stats group, there should be no change involving this.

Pending

• support storage of multiple observed values
  • cmdstan
  • cmdstanpy
  • numpyro
  • pyro
• update ic functions to use the data available
  • add a var_name/obs_name to waic, loo, compare...
  • take into account diferent cases such as an ic value per observation, considering all observations a single variable...
• add new group to inference data schema

This pull request adds a user guide for showing current dask capabilities in ArviZ. Closes #1631

Link to the notebook

• [x] Follows official PR format
• [ ] Changes are listed in changelog

Edits to be made:

• [x] Note on out of memory computation
• [x] Provide a link to best practices doc.
• [x] Section on how to get dask backed InferenceData objects.
• [x] For both non-dask and dask enabled method, initialize the data at the start instead of loading it everytime.
• [x] Small note on Client and ResourceProfiler
• [x] High level explanation on az.Dask and link to it's api docs
• [x] Check if dask="allowed" works or not

Hello, I am trying to load a pretty big (1.1GB) cmdstanpy csv output file using arviz with: data = az.from_cmdstan(posterior="path-to-csv") or with data = az.from_cmdstan(posterior=["path-to-csv", "parameter-name"]) both ways take more than an hour and use many GB of memory (at this point even 20gb is not enough and python process crushes). Is there a more efficient way to load the sampling output data?

This issue was also discussed here - https://discourse.mc-stan.org/t/loading-cmstanpy-output-in-python/17828

arviz - the current master version (0.9.0) cmdstanpy - 0.9.63 cmdstan - 2.24.1 python - 3.7.0 centOS7

Additional minimal documentation corrections for:

1 auto_corr 2 bpv 3 ess

plots

Checklist

• [x] Follows official PR format
• [x] Code style correct (follows pylint and black guidelines)
• [ ] Changes are listed in changelog

For feature parity with 1D KDE, custom limits could be supported in 2D.

I hacked an implementation by simply adding optional xmin, xmax, ymin, ymax arguments to _fast_kde_2d to override the limits inferred from the input data. Here's an example (left is current behavior of plot_kde, right has the same grey contours but the filled contour plot modifies xmin and ymin).

Just thought I'd throw this out there since the implementation is simple. Happy to share a reproducer (or draft a PR) if desired.

A kde_kwargs parameter to plot_kde, passed on to the KDE routines, would be useful. It seems the 2D KDE method is undocumented/private, though - making it public would be nice rather than restricting the flag to 1D KDEs only (though _fast_kde_2D has few parameters for customization). (IMHO, It would even nicer if kde could mimic plot_kde's generality to 1D and 2D by accepting a values2 argument, but that would be a breaking change since kde currently accepts circular as the second positional argument.)

The docstring of kde cites kde_linear and kde_circular for details on keyword args, but their docstrings don't appear in the documentation. (Note that both methods actually have leading underscores.)

Short Description

Greetings all,

I'm looking for a way to plot differences based on random effects on categorical regression. However, I failed to find any examples or API reference.

For instance;

Code Example or link

import arviz as az
import bambi as bmb
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
import pymc as pm
import pymc.sampling_jax
import numpyro

df = pd.DataFrame({"Target": target_variable, "Random_Effects_Var": random_effects, "Column1": column1, "Column2": column2})

model = bmb.Model("Target~ Column1 + Column2 + (1|Random_Effects_Var)",
                  data = df, family = "categorical", auto_scale = True, priors = priors)

model.build()

with model.backend.model:
    idata=pm.sampling_jax.sample_numpyro_nuts(draws= 250, tune = 50, target_accept = .99,
                                              postprocessing_backend = 'cpu', idata_kwargs={"log_likelihood":True},chain_method='parallel')
    posterior_predictive = pm.sample_posterior_predictive(trace = idata, extend_inferencedata=True)

idata

So far, so good. Imagine, Random_Effect_Var includes 2 levels of categorical data, A and B , and I'd like to plot the influence of each (Indeed group-level scores) over the Target variable along with variables in Column1 and Column2.

here is an example from R package Modelbasedlink

How can I plot something similar with Arviz?

Description

There was an issue with scipy 1.10 rc1, so I added this ravel to get CI green quickly. It should be fixed already in rc2 so there is no need to keep it.

Checklist

• [x] Does the PR follow official PR format?
• [x] Is the code style correct (follows pylint and black guidelines)?

:books: Documentation preview :books:: https://arviz--2187.org.readthedocs.build/en/2187/