Recent content on Dr. Juan Camilo Orduz
This example is part of the NumPyro Forecast documentation; see here . 
 This notebook extends that example by training on the full dataset ( \(50K\) time series) on a GPU via Modal .
The end-to end notebook runs in approximately \(10\) minutes on GPU 🚀! 
 
 The FreshRetailNet-50K dataset ( Yang et al., 2025 ) contains \(50{,}000\) daily store-product demand series from fresh…
In this notebook, we work out an example of how to interpret and communicate statistical models. We follow the ideas and techniques from the amazing book “Model to Meaning: How to Interpret Statistical Models with marginaleffects for R and Python” . This exposition is by no means exhaustive, but it should give you a good starting point. For more details, check the book! 
 
 Motivating…
In the blog post “Notes on Exponential Smoothing with NumPyro” we implemented the damped Holt-Winters model in the component form of Hyndman & Athanasopoulos, FPP3 , §8.3 . A very nice observation was raised by @theorashid in issue #193 :
 when we forecast with the component-form model the credible intervals barely widen across the horizon . The reason is subtle but important. Quoting the…
In this notebook we provide an introduction to the frugal parameterization for Bayesian causal inference, developed in the paper “Parameterizing and simulating from causal models”, (Evans & Didelez, 2024) . We illustrate the concepts with a simple example using synthetic data and an explicit PyMC implementation. 
 The authors of the frugal parameterization (Evans & Didelez, 2024) have…
Mediation analysis goes beyond asking “does the treatment work?” to ask “how does the treatment work?” Understanding the mechanisms by which an intervention achieves its effect can have important implications for what treatments or policy changes are preferable. For instance, a family intervention program during adolescence might reduce substance use disorder in young adulthood, but through which…
In many forecasting applications, probabilistic methods have become a standard tool. The reason is that predicting the point forecast is not enough for decision-making. The predicted uncertainty provides valuable information that can be used to make better decisions, such as risk management or inventory control. How to evaluate the performance of probabilistic forecasts? A naive approach is to…
This notebook implements a simulation-based Bayesian power analysis for conversion-rate A/B tests and compares it against the classical frequentist
two-proportion z-test. We follow the Kruschke HDI+ROPE framework from The Bayesian New Statistics (Section “Planning for precision and other goals: Better done Bayesian”), in the same way as in the introductory example: Introduction to Bayesian…
In many years of professional data science, I have seen many senior decision makers choose models based on in-sample \(R^2\) — the more the better. Many of them would not even tolerate a model with an \(R^2\) lower than \(0.8\) . These decision makers base their decisions on hyper inflated \(R^2\) models of the orders of millions of dollars (think about marketing spend). Here is a hard pill to…
In this notebook we reproduce the fantastic material from the video Statistical Rethinking 2026 Lecture B04 - Group-level confounding and intro to social networks by Richard McElreath. This is a great video to understand the difference between fixed and random effects models. It is a must watch! Here we use PyMC to fit the models and compare the results. 
 
 
 
 
 
 The…
In this notebook I explore the Bayesian Decision Theory Workflow described in the amazing blog post A Bayesian Decision Theory Workflow by Daniel Saunders . It is a great resource to understand how to use Bayesian methods for optimization problems. In the original post, Daniel explains the theory and key concepts using PyMC and its computational backend PyTensor . He does a remarkable job…
In this third notebook, I extend the hierarchical forecasting model from Part II by adding a neural network component to the state transition function. This creates a Hybrid Deep State-Space Model that combines probabilistic modeling with deep learning. 
 Why? This is a personal experiment to explore how to integrate neural networks with hierarchical models. It is not adding complexity for the…
Based on my experience, whenever someone asks for a prediction (or forecasting) model, they actually do not need a prediction model per se. They typically want to answer causal questions or do some kind of optimization. I have covered some case studies on causal questions in previous posts (for example, Introduction to Causal Inference with PPLs and “Using Data Science for Bad Decision-Making: A…
In this notebook, we demonstrate how to estimate Conditional Average Treatment Effects
(CATE) using a Causal Effect Variational Autoencoder (CEVAE) by implementing an example
from scratch in NumPyro . This approach is particularly useful when we suspect
the presence of unobserved confounders that affect both treatment assignment and outcomes. 
 Disclaimer : I am not an expert in…
This notebook demonstrates how to estimate the Average Treatment Effect (ATE) using variational
inference in the presence of unobserved confounders. The approach is based on the tutorial by Robert Ness from his book: Causal AI book . We port his Pyro code to NumPyro (this was the objective for me to learn the details of the method). 
 The fundamental challenge in causal inference from…
Estimating causal effects from clustered or grouped data requires careful attention to the hierarchical structure of observations. When units are nested within groups such as students within classrooms, or patients within hospitals—ignoring this structure can lead to incorrect standard errors, inefficient estimates, and invalid causal inferences. Multilevel models provide a principled framework…
Causal inference asks a deceptively simple question: “What would have happened if things were
different?” Whether we’re evaluating a job training program, testing a new medical treatment,
or analyzing the impact of a policy change, we want to understand the causal effect of an
intervention not just observe correlations in the data. 
 Traditional statistical methods often struggle…
In this notebook, we present how to implement and fit Bayesian Vector Autoregressive (VAR) models using NumPyro . We work out three components: 
 
 Specifying and fitting the model in NumPyro 
 Using the model to generate forecasts 
 Computing the Impulse Response Functions (IRFs) 
 
 We compare these three components with the ones obtained using the statsmodels…
In this notebook we provide a brief introduction to Stochastic Variational Inference (SVI) with NumPyro . We provide the key mathematical concepts, but we focus on the code implementation. This introductory notebook is meant for practitioners. We do this by working through two examples: a very simple parameter recovery model and a Bayesian Neural Network. 
 This work was presented at PyData…
In this notebook, we extend the revenue-retention model introduced in the sequence of blog posts “Cohort Revenue & Retention Analysis: A Bayesian Approach” and “Cohort Revenue Retention Analysis with Flax and NumPyro” (plus the associated pre-print “Cohort Revenue & Retention Analysis: A Bayesian Approach” ) to include analysis across different markets (or any type of grouping variable). The…
In this notebook, we present a simple example to illustrate how to vectorize the ROC curve computation over a Bayesian model. This is helpful when we want to compute the ROC curve using the implementations from scikit-learn . We use the classical moons dataset to generate the data and fit a Gaussian process models similarly as the previous post Scikit-Learn Example in PyMC: Gaussian Process…
Recently, I have been thinking a lot about data-driven decision-making, particularly in the context of experimentation. Why? I am uncomfortable with the common practice of using p-values and frequentist null hypothesis significance testing to make decisions. I don’t feel confident about the approach. I think it is because I do not get it. For instance, when I am forced to explain the definition of…
In this notebook we generate a simulation to reproduce the results of the great blog post “The Bet Test: Spotting Problems in Bayesian A/B Test Analysis” , where Tyler Buffington discussed about some caveats of using Bayesian AB testing when not thinking about the prior predictive distribution on the key metrics of interested. If you haven’t read it yet, please do it before continuing ;) 
 In…
In this notebook, we want to explore some ideas on hierarchical Hilbert Space Gaussian Processes following the fantastic exposition of the PyMC example notebook: “Gaussian Processes: HSGP Advanced Usage” by Bill Engels , Alexandre Andorra and Maxim Kochurov . I can only recommend to read the notebook and the references therein!
For an introduction to Hilbert Space Gaussian Processes, please…
In many demand forecasting problems, we face the challenge of predicting the demand for a product or service with very sparse data. This is especially true in retail, where many products are only sold occasionally. This sparsity can come from many sources, as greatly described in the (recommended!) blog post “Why zeroes happen” by Ivan Svetunkov . On one hand, the sparsity can simply come from the…
Motivated by the great blog post by Understanding CUPED by Matteo Courthoud , we explore a Bayesian approach to CUPED to understand its sensitivity with respect to the covariance parameter \(\theta\) of the pre-post mean model. I will assume that the reader is already familiar with CUPED and has read Matteo’s blog post (highly recommended!). Here we focus on the sensitivity component. We do not do…
We present an example of a dynamic forecasting time-series model that incorporates a prior calibration process to estimate the temperature effect on electricity demand. The model is based on the previous example Electricity Demand Forecast: Dynamic Time-Series Model . In this second iteration, we borrow the ideas from the Pyro great example Forecasting with Dynamic Linear Model (DLM) where they…
We work out a classical electricity demand forecasting model form the case study Structural Time Series Modeling Case Studies: Atmospheric CO2 and Electricity Demand from the TensorFlow Probability documentation. The idea of this example is to use temperature as a linear covariate to model the electricity demand. In this example, we show how to use a (Hilbert Space Approximation) Gaussian process…
In this second notebook, we continue working on the NumPyro implementation of the hierarchical forecasting models presented in Pyro’s forecasting documentation: Forecasting III: hierarchical models . In this second part, we extend the model described in the first part From Pyro to NumPyro: Forecasting Hierarchical Models - Part I by adding all stations to the model. 
 
 Prepare Notebook…
In this notebook we provide a NumPyro implementation of the first model presented in the Pyro forecasting documentation: Forecasting III: hierarchical models . This model generalizes the local level model with seasonality presented in the univariate example Forecasting I: univariate, heavy tailed (see From Pyro to NumPyro: Forecasting a univariate, heavy tailed time series for the corresponding…
In this notebooks we port the Pyro forecasting example Forecasting I: univariate, heavy tailed to NumPyro . The forecasting module in Pyro is fantastic as it provides an easy interface to develop custom forecasting models. It has also many helpful utility functions to generate features and for model evaluation. The purpose of this translation is to dig deeper into the some forecasting components…
In this notebook we use a retail publicly available dataset to fit and compare various pricing elasticity models. This example can be seen as a continuation of the notebooks regarding Bayesian hierarchical models (see for example Multilevel Elasticities for a Single SKU - Part I ). In this example we also see how Bayesian hierarchical models can help regularize elasticity estimates when taking…
In this notebook we continue our simulation study for elasticities, see: 
 
 Multilevel Elasticities for a Single SKU - Part I 
 Multilevel Elasticities for a Single SKU - Part II 
 
 for an introduction to the problem and some models. We extend the covariance model to allow two covariance components on both the intercepts and slopes (coefficient of the median_ income…
In this blog post, we experiment with a hierarchical exponential smoothing forecasting model, extending the ideas from the univariate case presented in the blog post “Notes on Exponential Smoothing with NumPyro” . We use NumPyro and compare the NUTS and SVI results. For such a purpose, we use Continuous Ranked Probability Score ( CRPS ). We also compare these forecasts with univariate statistical…
In this notebook we will explore the use of censored likelihoods for demand forecasting. 
 Business Problem : Let us assume we have a store with a single product. We want to forecast the (true!) demand for this product for the next 30 days using historical data. The historical data consists of the daily sales of the product for the last year (approximately). An important challenge is that our…
In this notebook, we explore the (conceptual) ideas and (practical) implementation details of the Hilbert Space approximation for Gaussian processes introduced in the article “Hilbert space methods for reduced-rank Gaussian process regression” by Arno Solin and Simo Särkkä. We do not go deep into the mathematical details (proofs) but focus on the core ideas to help us understand the main concepts…
In this notebook, we explore how we can use Bayesian modeling to estimate the parameters of a censored data set. These datasets are common in many fields, including survival analysis and supply chain management. I was motivated to write this notebook after reading the excellent blog post “Modeling Anything With First Principles: Demand under extreme stockouts” by Kyle Caron where he uses these…
After going through the fundamentals of the TSB Method for Intermittent Time Series Forecasting in NumPyro in the previous notebook, we explore a variation of it that might be useful for certain applications. In a nutshell, we keep the same model structure of the TSB model, but we modify the likelihood function to account for the sparsity of the time series. Concretely, we replace the classic…
In this notebook we provide a NumPyro implementation of the TSB (Teunter, Syntetos and Babai) method for forecasting intermittent time series. The TSB method is similar to the Croston’s method in the sense that is constructs two different time series out of the original one and then forecast each of them separately, so that the final forecast is generated by combining the forecasts of the two time…
In this notebook, we will implement Croston’s method for intermittent demand forecasting using NumPyro . Croston’s method is a popular forecasting method for intermittent demand data, which is characterized by a large number of zero values. The method is based on the idea of separating the demand size and the demand interval, and then forecasting them separately using simple exponential smoothing.…
This are some notes on how to implement an ARMA(1, 1) model using NumPyro for time series forecasting. The ARMA(1, 1) model is given by 
 \[y_t = \mu + \phi y_{t-1} + \theta \varepsilon_{t-1} + \varepsilon_t\] 
 where \(y_t\) is the time series, \(\mu\) is the mean, \(\phi\) is the autoregressive parameter, \(\theta\) is the moving average parameter, and \(\varepsilon_t\) is a white noise…
This notebook serves as personal notes on NumPyro ’s implementation of the classic exponential smoothing forecasting method. I use Example: Holt-Winters Exponential Smoothing . The strategy is to go into the nitty-gritty details of the code presented in the example from the documentation: “Example: Holt-Winters Exponential Smoothing” . In particular, I want to understand the auto-regressive…
In this notebook, we present a complete simulation study of the media mix model (MMM) and experimental calibration method presented in the paper “Media Mix Model Calibration With Bayesian Priors”, by Zhang, et al. , where the authors propose a convenient parametrization the regression model in terms of the ROAs (return on advertising spend) instead of the classical regression (beta) coefficients.…
In this notebook we present an alternative implementation of the cohort-revenue-retention model presented in the blog post Cohort Revenue & Retention Analysis: A Bayesian Approach where we show how to replace the BART retention component with a general neural network implemented with Flax . This allows faster inference, as we can use NumPyro ’s NUTS sampler or any of the stochastic variational…
In this notebook I want to experiment with the numpyro/contrib/module.py module which allow us to integrate Flax models with NumPyro models. I am interested in this because I want to experiment with complex bayesian models with larger datasets. 
 Most of the main components can be found in the great blog post Bayesian Neural Networks with Flax and Numpyro . The author takes a different path…
In this notebook we want to reproduce a classical example of using Gaussian processes to model time series data: The birthdays data set. I first encountered this example in the seminal book Chapter 21, Bayesian Data Analysis (Third edition) when learning about the subject. One thing I rapidly realized was that fitting these types of models in practice is very computationally expensive and…
In this notebook we present an example of how to use a combination of Bayesian hierarchical models and the non-parametric methods , namely bayesian additive trees ( BART ), to model the product life cycles. This approach is motivated by previous work in cohort analysis, see here . 
 As a case study we use the Google search index (trends) data for iPhones worldwide. We use the data of four…
In this notebook we describe how to use blackjax ’s pathfinder implementation to do inference with a numpyro model. 
 I am simply putting some pieces together from the following resources (strongly recommended to read): 
 References: 
 
 Blackjax docs: Use with Numpyro models 
 Blackjax Sampling Book: Pathfinder 
 Numpyro Issue #1485 
 PyMC Experimental - Pathfinder…
I had the great opportunity to have a conversation with Aleksander Molak about causality, marketing and simulations as well as my career from academia to industry. Check it out! 
 
 
 

In this notebook we go deeper into the last covariance model presented in the previous blog post Multilevel Elasticities for a Single SKU . In particular we describe how to generate posterior predictive samples from an unseen region by the model. This can be useful for scenario planning: once can simulated outcome quantities from price ranges through the elasticity estimates (with uncertainty!)…
In this notebook I want to experiment with some basic models for price elasticity estimation in the simple context of a simple sku across multiple stores and regions. The motivation is to have a concrete example of the elasticity models presented in the Chapter 11: Big Data Pricing Models of the book Pricing Analytics by Walter R. Paczkowski . 
 
 Elasticity Definition 
 Here I provide…