Proudly sponsored by PyMC Labs, the Bayesian Consultancy. Book a call, or get in touch!

• My Intuitive Bayes Online Courses
• 1:1 Mentorship with me

What’s the difference between MCMC and Variational Inference (VI)? Why is MCMC called an approximate method? When should we use VI instead of MCMC?

These are some of the captivating (and practical) questions we’ll tackle in this episode. I had the chance to interview Charles Margossian, a research fellow in computational mathematics at the Flatiron Institute, and a core developer of the Stan software.

Charles was born and raised in Paris, and then moved to the US to pursue a bachelor’s degree in physics at Yale university. After graduating, he worked for two years in biotech, and went on to do a PhD in statistics at Columbia University with someone named… Andrew Gelman — you may have heard of him.

Charles is also specialized in pharmacometrics and epidemiology, so we also talked about some practical applications of Bayesian methods and algorithms in these fascinating fields.

Oh, and Charles’ life doesn’t only revolve around computers: he practices ballroom dancing and pickup soccer, and used to do improvised musical comedy!

Our theme music is « Good Bayesian », by Baba Brinkman (feat MC Lars and Mega Ran). Check out his awesome work at https://bababrinkman.com/ !

Thank you to my Patrons for making this episode possible!

Yusuke Saito, Avi Bryant, Ero Carrera, Giuliano Cruz, Tim Gasser, James Wade, Tradd Salvo, William Benton, James Ahloy, Robin Taylor,, Chad Scherrer, Zwelithini Tunyiswa, Bertrand Wilden, James Thompson, Stephen Oates, Gian Luca Di Tanna, Jack Wells, Matthew Maldonado, Ian Costley, Ally Salim, Larry Gill, Ian Moran, Paul Oreto, Colin Caprani, Colin Carroll, Nathaniel Burbank, Michael Osthege, Rémi Louf, Clive Edelsten, Henri Wallen, Hugo Botha, Vinh Nguyen, Raul Maldonado, Marcin Elantkowski, Adam C. Smith, Will Kurt, Andrew Moskowitz, Hector Munoz, Marco Gorelli, Simon Kessell, Bradley Rode, Patrick Kelley, Rick Anderson, Casper de Bruin, Philippe Labonde, Michael Hankin, Cameron Smith, Tomáš Frýda, Ryan Wesslen, Andreas Netti, Riley King, Yoshiyuki Hamajima, Sven De Maeyer, Michael DeCrescenzo, Fergal M, Mason Yahr, Naoya Kanai, Steven Rowland, Aubrey Clayton, Jeannine Sue, Omri Har Shemesh, Scott Anthony Robson, Robert Yolken, Or Duek, Pavel Dusek, Paul Cox, Trey Causey, Andreas Kröpelin, Raphaël R, Nicolas Rode, Gabriel Stechschulte, Arkady, Kurt TeKolste, Gergely Juhasz, Marcus Nölke, Maggi Mackintosh, Grant Pezzolesi, Avram Aelony, Joshua Meehl, Javier Sabio, Kristian Higgins, Alex Jones, Gregorio Aguilar and Matt Rosinski.

Visit https://www.patreon.com/learnbayesstats to unlock exclusive Bayesian swag ;)

Links from the show:

• Charles’ website: https://charlesm93.github.io/
• Charles on Twitter: https://twitter.com/charlesm993
• Charles on GitHub: https://github.com/charlesm93
• Charles on Google Scholar: https://scholar.google.com/citations?user=nPtLsvIAAAAJ&hl=en
• Stan software: https://mc-stan.org/
• Torsten – Applications of Stan in Pharmacometrics: https://github.com/metrumresearchgroup/Torsten
• R̂ – Assessing the convergence of Markov chain Monte Carlo when running many short chains: https://arxiv.org/abs/2110.13017
• Revisiting the Gelman-Rubin Diagnostic: https://arxiv.org/abs/1812.09384
• An importance sampling approach for reliable and efficient inference in Bayesian ordinary differential equation models: https://arxiv.org/abs/2205.09059
• Pathfinder – Parallel quasi-Newton variational inference: https://arxiv.org/pdf/2108.03782.pdf
• Bayesian workflow for disease transmission modeling in Stan: https://mc-stan.org/users/documentation/case-studies/boarding_school_case_study.html
• LBS #76 – The Past, Present & Future of Stan, with Bob Carpenter: https://learnbayesstats.com/episode/76-past-present-future-of-stan-bob-carpenter/
• LBS #51 – Bernoulli’s Fallacy & the Crisis of Modern Science, with Aubrey Clayton: https://learnbayesstats.com/episode/51-bernoullis-fallacy-crisis-modern-science-aubrey-clayton/
• Flatiron Institute: https://www.simonsfoundation.org/flatiron/
• Simons Foundation: https://www.simonsfoundation.org/

Abstract

by Christoph Bamberg

In episode 90 we cover both methodological advances and their application, namely variational inference and MCMC sampling and their application in pharmacometrics. 

And we have just the right guest for this topic - Charles Margossian! You might know Charles from his work on STAN, his workshop teaching or his work at his current position at the Flatiron Institute.

His main focus now is on two topics: variational inference and MCMC sampling. When is variational inference (or approximate Bayesian methods) appropriate? And when does it fail? Charles answers these questions convincingly, clearing up some discussion around this topic.

In his work on MCMC, he tries to answer some fundamental questions: How much computational power should we invest? When is MCMC sampling more appropriate than approximate Bayesian methods? The short answer: when you care about quantifying uncertainty. We even talk about what the R-hat measure means and how to improve on it with nested R-hats.

After covering these two topics, we move to his practical work: pharmacometrics. For example, he worked on modelling the speed of drugs dissolving in the body or the role of genetics in the workings of drugs. 

Charles also contributes to making Bayesian methods more accessible for pharmacologists: He co-developed the Torsten library for Stan that facilitates Bayesian analysis with pharmacometric data. 

We discuss the nature of pharmacometric data and how it is usually modelled with Ordinary Differential Equations. 

In the end we briefly cover one practical example of pharmacometric modelling: the Covid-19 pandemic.

All in all, episode 90 is another detailed one, covering many state-of-the-art techniques and their application.

Transcript

This is an automatic transcript and may therefore contain errors. Please get in touch if you're willing to correct them.

Charles Margossian, welcome to Learning Vision Statistics.

Hi Alex, thanks for having me.

Yeah, thanks a lot for taking the time. I'm super happy to have you here.

It's been a while since I wanted to have you on the show and now I managed to

find the slot and thank you also to a few patrons who sent me messages

to tell me that they would really like to hear about you on the show. So thanks

a lot folks for being so proactive and giving me. ideas for the show. So Charles,

let's dive in and

But

as

maybe

usual...

I can bring up a very quick anecdote, which

Oh sure.

is I think two, three weeks ago, your show came up. And I told my

colleagues, I would feel like I have made it in the Bayesian. statistics

world if I get an invitation to speak on learning Bayesian statistics. So

I'm thrilled to be here and I think sort of the mysterious patrons

who have incited this meeting.

Yes, I'm sure they will recognize themselves. So, yeah, let's start with your

origin story, actually. That's also something I found super interesting. How

did you come to the world of statistics and pharmacometrics and epidemiology?

Mm-hmm. Right. So I became interested in statistics and data as an undergrad.

Actually,

Mm-hmm.

I was studying physics. I was already in the US. I was at Yale and

I was working with an astronomy lab

Mm-hmm.

on Exoplanet. We had a lot of data.

Mm-hmm.

So that generally got me interested in data science. But I wasn't

introduced to Bayesian methods. And I always was a bit uncomfortable

with how we handling uncertainty, what might we do about it. And I was fortunate

after I graduated to actually get a job at a, at a biotech company where,

so the company was met from research group. They were

Mm-hmm.

based in Connecticut and the supervisor I got there, William Gillespie.

Mm-hmm.

is one of the pioneers of Bayesian methods in pharmacometrics. And

so he introduced me to Bayesian statistics and it made a lot of sense.

And I realized, oh, this is what I have been looking for when I was

doing astronomy. Right. And I started, you know, piecing this together.

And to be fair, I hadn't done a lot of statistics before. I don't have

the experience of struggling with classical statistics for a decade before

being rescued by Bayesian statistics. I encountered Bayesian very early

on. But even then, the way we counterfeit uncertainty, the fact that

everything reduced to essentially one equation was extremely compelling

to me. And Bill Gillespie, hired me

Mm-hmm.

at Matrim Research Group to work on Stan. So this was back in 2015.

And at the time, Stan was fairly new, very exciting, very promising.

And you had these specialized software in pharmacometrics, but they were

not open source. They were not really supporting Bayesian statistics,

or at least that was not their priority. And they didn't give you flexibility

that a lot of my colleagues were looking for. On the other hand, you had

Stan, which was open source, which had great algorithms to do vision modeling,

but which lacked a lot of the features required to do pharmacometrics

modeling, right? One example was support for differential equation solvers.

And all these models are based

Yeah.

on ODEs. And at the time, Stan had limited support for ODEs. and more

generally implicit functions.

And then there were some more specialized things like handling, you

know, the event schedule of clinical trials. And the project ended up

being, we're going to write some general features, we're going to

contribute them to Stan, and then we're going to write a specialized

extension called Torsten. And that's going to have, you know, more bespoke

features targeted at a pharmacometrics audience. And so... That was

my exposure to Stan. I can tell you about my first poll request, which

Bob Carpenter and Daniel Lee reviewed

Mm-hmm.

extensively. They were very patient with my C++. And

Hehehe

my plan, my loose plan had been, well, I'm gonna work for one or two

years. I'm gonna learn programming and statistics because I had gained

an appreciation for it in astronomy. And then I'm going to do a PhD

in physics. But Andrew Gellman, who I had met in an official capacity

like through my work, encouraged me to apply to the statistics program

at Columbia University, saying, you know, if you do statistics, you'll

still be able to work on the natural sciences and the physics that

you're interested in, but you

Yeah.

know, you'll have, maybe you'll be able to make a more unique contribution

as a statistician. So I took his word for it. And yeah, and once

I had the offer from Columbia, I mean, it was a difficult decision. Um, but

I decided, okay, I'm going to do a PhD in statistics. And that was a big

change, the field. Um, but I did pursue that. I was able to continue

working on Stan and on Stan adjacent, uh, projects. And just a year

ago, I completed the PhD. So I would say that's the origin story. Um,

that eventually led me to the Flatiron Institute. So that's where

I'm currently, it's in New York. We're a non-profit. We focus on applying

computational methods to the basic sciences.

Mm-hmm.

So big emphasis on collaboration. We have people who do astronomy, who do quantum

physics, do biology. And so that really resonates with why I wanted

to do statistics in the first place, which I wanna solve problems in

the sciences.

Yeah, I can really relate to that. That was also kind of what happened to

me too, even though the first field was political science and not astronomy.

Okay.

But yeah,

Yeah.

definitely in the end, applying the methods became more interesting than the

field in itself. So continued

Right,

on that path.

right, and I think that, you know, the attitude I had, so I had offers

to do a PhD in physics, a PhD in statistics, and a bunch of different

fields. And at the time, my attitude was, look, wherever I go, I'm gonna

have to do, you know, analyze some data, understand where that data

comes from. I'm gonna have to, you know, have the tools to do the

analysis or do the statistics.

Mm-hmm.

So the idea is, I think our loyalty is not to a discipline or to a

field, it's really to a problem.

Mm-hmm.

Whatever

Yeah.

problem we work on, there's always a lot that we need to learn. We're never

experts in a new problem, in a new research problem. And so kind of like,

try not to take the field that I mean, whether it's for the PhD or even,

right now I'm in computational mathematics, right? You know, whatever that

means. that I care more about the problems I work on than the department

I'm affiliated with.

Yeah, yeah, completely, completely really do that for sure. I like the fact of

using the method to solve an interesting problem, whether it is about astronomy,

political science, biology,

Mm-hmm.

epidemiology, it's like, and in the end, I find that even more interesting

because you get to work on a variety of topics that you wouldn't have otherwise.

And also you

Right.

learn so much. So that's really

Yeah,

the cool thing.

yeah.

Yeah. Thanks a lot for this introduction. And so right now, well, these

days you're working at the Flatiron Institute, as you were saying, and it seems

to be a very diverse organization. We had Bob Carpenter actually

Mm-hmm.

recently on the show, so I'm going to link to his episode in the show notes.

But you, Charles, what are the topics you are... particularly interested in

these days at the flat iron.

Yeah, I think that I have a little bit two pools right now. There are

a bit on the methodology side

Mm-hmm.

of things. And one of them is variational inference, which is a form

of, it's called approximate Bayesian method, is what it's called. And really

what it's... What I'm trying to understand is when should we use

an approximate method like variational inference, because it's incredibly

popular. It's used in a lot of fields of machine learning.

Mm-hmm.

And it's used a lot of times in artificial intelligence. And yet you

have so many smart, brilliant people who are completely distrustful of variation

and inference. And it's also not difficult to construct an example

where it really doesn't do the job, where it really fails in a spectacular

way. And I think that ultimately it depends on the problem you apply

it to. And what we need to do is understand. you know, when can we

get away with the approximations that variational inference proposes to do?

Why does it sometimes work really well? Or why do we sometimes really

get punished for using variational inference? And I can give you a

very simple example of that,

Yeah.

which was, so this was a recent work with. with Lauren Saul, who is

also at the Flatiron Institute. He kind of has machine learning there. And

there's a statement that goes around, which is to say, well, variation

inference will give you good estimates of the expectation value for

a target distribution. So that's great. We often care about expectation

values, whether

Mm-hmm.

that's in statistical physics or in Bayesian statistics.

Mm-hmm.

But on the other hand, it will underestimate uncertainty.

Mm-hmm. Okay.

Okay, well, what does it mean to underestimate uncertainty? How do

we, you know, how do we make sense of a statement like that, which

has appeared time and time again, over the last two decades.

Mm-hmm.

And really we realized that there were two measures of uncertainty

that seem to come up again and again, one is the marginal variances and

the other one is the entropy

Mm-hmm.

and so people like the entropy because it's a multivariate notion.

of uncertainty and it's a generalization of the marginal variance.

So when people say variation on inference underestimates uncertainty, they

usually mean you're underestimating the marginal variances and you're underestimating

the entropy.

Okay.

And so what we ended up doing is demonstrating this on the Gaussian

case, right? You have a Gaussian target, you're approximating it

with a Gaussian with a diagonal covariance matrix. So that's called

the factorized approximation or the mean field approximation.

Mm-hmm.

And indeed you underestimate the marginal variance and you underestimate

the entropy. Here's where it gets interesting, which is as you start

going to higher dimensions, And for some reason, I don't think people,

people have only looked at, you know, really the two dimensional case,

because that's the figure that fits on that page. Right. But if you start

going to higher dimensions and you take limits where the dimension goes

to infinity, you can construct examples where you actually get very, very

accurate estimates of the entropy, but you are underestimating the

marginal variances in every dimension in an arbitrarily bad man.

And so the two notions of uncertainty are not at all equivalent and

not at all interchangeable. And what ends up happening is you look at fields

where variation inference is applied. So

Mm-hmm.

for example, in statistical physics, where, where people want to estimate

the entropy of an easing model, for example, and that's

Mm-hmm.

where this factor is this mean field approximation comes

Okay.

from. Well, here it works just fine. You actually get good estimates

of the entropy, right? In certain limits. In machine learning, where

you're trying to maximize marginal likelihoods, rather than default Bayesians,

actually, you know, you get good estimates of those marginal likelihoods.

At least that's our working conjectures. But in

Okay.

Bayesian statistics, where we have interpretable quantities, and we know

those marginal variances mean something for those parameters that have

a meaning. Well, here, variation inference might really not or at

least not vanilla implementations in it. And that's an example of how, by studying

an example, we can start understanding why is it that some people are

so enthusiastic about variation inference and other people are so

distrustful of it. It really depends on what is the measure of uncertainty

that you care about, and that in turn is informed by the problem you

want to solve. So I bring this up as an archetype of the work that we're

trying to do.

Mm-hmm.

We wanna understand this method one in one.

Mmm. Yeah, that's really interesting.

First big topic.

Yeah.

The other big topic is I still do a lot of MCMC.

Mm-hmm.

I care a lot about MCMC because as we'll see in pharmacometrics, I

really think that is our best tool to solve many problems.

Mm-hmm.

And here I'm trying to understand some more fundamental questions,

right? So people often say, well, you know, yeah, MCMC is great but

it's too computationally expensive. That's why we'll

Yeah.

use approximate methods. And one thing I will argue is actually it's

computationally expensive, but we don't really have a very good sense

of how much computation we should throw at MCMC. Because ultimately

we have three fundamental tuning parameters. One is the number of chains

that we use. The other one is how long is the warmup or the burning

phase? And then the third one is how long is the sampling phase? And

actually,

it's not clear what is the optimal computation that you should throw

at an MCMC problem. I think people rely on heuristics. More often,

they rely on conventions. These are the defaults in Stan, in PMC, in

TensorFlow probability. But actually, you know, we need to think about

have we used a warm up phase that's too long or too short or

Mm-hmm.

sampling phase that's too long or too short, right?

Yeah.

Or how many chains is it useful to do? Especially now we have GPUs,

right? In theory, I could run MCMC with 10,000 chains. And people do

that now, right? I mean, some people do that. It selects you do that

now, right? But what are the implications of doing that, right? I'll

tell you an implication. If I have 10,000 chains, that I'm running on

a GPU. then I think the sampling phase can be one iteration per chain.

We can discuss this and we can argue for or against it, but it changes

a little bit our perspective of how much computation we throw at it. Another

perspective I like is if you have a strict computational budget,

and then you say, well, I'm not going to use MCMC, it's too expensive.

Let me run variation inference. I know it's biased. I know it's approximate,

but at least it, it finishes running in, you know, within 10 minutes and

like, well, run 10 minutes of MCMC. It'll be biased. It'll be approximate

and it will finish running in 10 minutes. And then ask yourself, well,

how good is this estimate? thinking a little bit more carefully about,

you know, how much computation do we really need for MCMC and for the different

problems that we might be trying to solve.

Mm-hmm.

Yeah, that's really interesting. Thanks a lot for that very clear presentation.

I really love it. It's actually continue on that path because I wanted to ask

you about that a bit later in the show anyways. Yeah, several things that

bumped into my mind. First thing is MCMC is also an approximation method. So...

Why

Yes.

do we say, why do we, and I know we usually do that in the field, we define

variational inference as an approximation method, which kind of underlies,

assume that MCMC is not, but it is. So can you maybe draw the distinction?

What makes the difference between the two methods? And why do we call them approximation

for variation?

Yeah, absolutely. So I think that. people think like a lot of statisticians

asymptotically. So you know that asymptotically in what sense when

you run a single chain for an infinite number of iterations, MCMC

is not only gonna generate samples from stationary distribution, which

oftentimes is the posterior distribution, but also multicolor estimators

with an arbitrary precision. Your multicolor estimator will converge true

expectation value or the true variance or whatever it is you're trying

to estimate. Whereas with variational inference, here we have to be a

little bit careful because what does it mean for the asymptotic of variation

inference? So you might say, okay, I'm gonna run the optimization for

an infinite number of iterations. So let's assume that the optimizer

does converge. And actually until recently, this was not really shown.

There's a recent pre-print that I know Robert Gower and Justin Dumkey

have been working on where they actually showed that, yes, under certain

conditions, stochastic optimization will converge for variation inference.

Mm-hmm.

But then even if it doesn't converge, you say, well, let me think about

the approximation that minimizes my objective function. So oftentimes

the Colbert library divergence. That's what I get asymptotically. Well, that

will still not be in general, my target distribution.

Mm-hmm.

Even asymptotically, I'm still approximate. I think that's why people

draw this distinction between MCMC and variational inference. It's

really in the asymptotic sense that MCMC is an exact method, whereas

VI remains approximate. Even if you've thrown an infinite amount of

computation, you're probably not. Right? Now, in practice, we're not asymptotic,

right? We work with finite computation. And so I think it's very important

to recognize that, yes, MCMC is also an approximate method. It's not

unbiased because you don't initialize from the stationary distribution.

You do not reach the stationary distribution, right? So when I hear statements

like, well, first we wait for MCMC to converge to the stationary distribution.

Mm-hmm.

I think people have the right intuition, right? I don't think it's

a misleading statement, but nonetheless, it's an incorrect statement.

What we do is we wait for MCMC to get close enough to their stationary

distribution. And the reason why we care about being close enough

to the stationary distribution is because we want the bias to be small, right?

And so

Thank you.

when we think about, you know, convergence diagnostic. the way I think

we really should start thinking about a quantity like the r hat statistics,

for example. And r hat is interesting. r hat has been around for three

decades and frankly, there's still debates about what does r hat measure.

I mean, this is very existential, right? We have an estimator, but it's not clear

what the estimate is. And my perspective, my most recent perspective

is that what really matters, the reason I care about convergence is because

I want the bias of my Monte Carlo estimator to be sufficiently small.

It's not going to be zero, but it has to be sufficiently small. And so

can R hat tell me something about how small my bias is? And here's the

paradox, which is that R hat the way you compute it, it's a ratio

of two standard deviations. Right. So you know that when you measure variance,

it doesn't tell you something about bias. Right. And yet, you know, we say

our hat tells you if your warmup phase is long enough.

Yep.

Not everyone agrees that that's what it tells you, but you know, that's

my perspective and, um, and I think it's a reasonable perspective,

right? But the, the point of the warmup phase, you know, at a, the primary

point, not the only point is for the bias to go down. Right? So we're

faced with this paradox, this very fundamental question. Can R hat actually

give us any useful information? And there was a recent paper that argued that,

well, since R hat is really just looking at the variance, it's a one-to-one

map. So this was some really nice work by,

I know that Kida Vats is one of the co-authors on the paper, and

there's another co-author. They call it revisiting the Gellman-Rubin statistic.

Mm-hmm.

And they argue, well, R hat is just a reframing of variance and of

effective sample size. And to me, effective sample size only matters

when you're looking at the sampling phase, because it tells you, is

your variance low enough? And so we had to think a little bit hard about

this because it wasn't completely satisfactory because either it means

that r hat is not a useful convergence diagnostic in some sense, or actually

there's more going on. And what we realize is you look at the variance

that's being measured by r hat. So really what r hat ends up measuring

is you're running a bunch of chains. It could be four chains, it could

be more. Each chain generates one Monte Carlo estimator. And then you average

the per chain Monte Carlo estimator. So now you look at the variance

of a single chain Monte Carlo estimator, that variance you can actually

decompose it into a non-stationary variance and a persistent variance. And what

we realize, and you have to be careful, you have to do that analysis

for non-stationary Markov chains. Otherwise

Mm-hmm.

you completely miss the non-stationary variance.

Mm-hmm.

And the non-stationary variance is, you know, actually, it's a measure

of how well you've forgotten your initial point.

Mm-hmm.

And you can show in some cases that it decays at the same rate as the

squared bias. So you're not directly measuring the squared bias, but

because NCMC is what it is, the non-stationary variance gives you a

proxy clock for the bias. And so our argument is that, well, what's interesting

about our hat is not that it measures, you know, the persistent variance,

which you can then relate to the effective sample size, but that it measures

the non-stationary variance, right? And this then led us to, you know, we're

coming up with revisions of our hat, which more directly measure the

non-stationary variance rather than the total variance. So that we

actually get an estimator that unambiguously tells you something about

the length of the warm up phase. And then you ask the question of

the length of the sampling phase in a second and separate state.

So, sorry, it's conceptual explanation. So that this is a paper that we

have a preprint out called the nested r hat.

Mm-hmm.

And this is joint work with, so Andrew Gellman, Akiva Terry, Matt

Hoffman, so some of the usual suspects. We also have Pavel Junssov

and Yonel Viyouduran.

Mm-hmm.

So we've all worked on this together. And

A cool anecdote here is that what we were really interested in is those

regimes where we're running hundreds of chains in parallel, or thousands

of chains in parallel, like GPU-friendly MCMC. And what this made

us do is instead of thinking asymptotics in the limit where we have an infinitely

long chains, which is how asymptotics for MCMC have worked for the

past, you know, five or six decades, right? because that's what people

did. They ran long chains, right?

Yep.

And so even though asymptotics are property of infinity, we want to somehow

get close to asymptotic regime, right? That's where we care about this

asymptotic analysis. And here we thought, well, let's take asymptotics

in another direction. Let's say we have a finite number of chains,

but what happens when we have an infinite number of chains? And then

suddenly you can do asymptotic analysis. on non-stationary Markov chains.

Right? So the problem is if I take an asymptotic in the length of

the chain, well, I've made my chain stationary. And then there are

only so many properties that I can study. If I take asymptotics in

the other direction, which is the number of chains, and

Mm-hmm.

suddenly I can start making statements about non-stationary Markov chains,

I can elicit terms such as, you know, non-stationary variance. And

this was a cool example where the hardware you know, trying to work

with, you know, GPUs running a lot of chains. Actually, I think is really

changing our theoretical approach and our conceptual understanding of

MCMC. And I think we're going to get a lot of breakthroughs from this

kind of perspective.

Hmm. Yeah, super interesting. I put the two papers you mentioned in the show

notes. And so that makes me think, basically, what would you say right

now practically for people, when would variational inference usually be

most helpful, especially in comparison to MCMC?

Yeah, so great. So I think there are two things to unpack here. One is

when can variational inference still give you accurate answers?

Mm-hmm.

And two, when do you actually not need an accurate answer?

Right? So, in the first question, it really depends on the family of

variation and inference that you're willing to use, the family of approximation.

And also it turns out the objective function. So I mentioned earlier

that when I have this mean field approximation, and I'll just remind

what mean field means is, I'm assuming that all my latent variables

are independent. Now we don't believe that that's true in practice,

but that makes the computation much cheaper. And when you have, millions of

observations, millions of parameters, you need an algorithm where

the cost scales linearly with the number of observations. And if I don't

have this mean field assumption, I get things that can scale quadratically or

cubically.

And so when I do this approximation, well, I'm gonna get some things

wrong, but maybe I'll get the things that I care about right, which

Mm-hmm.

could be the first moment or it could be the entropy, right? if you change

the objective function, which actually the KL divergence, you can kind

of reverse it. You can show that you get arbitrarily poor estimates

of the entropy, but good estimates of the marginal variances. Right,

so it turns out that the choice of objective function that you use to

measure, you know, the disagreement between your approximation and

your target matters. So there's a real question of, you know, what

are the quantities you care about? because we don't care about the

whole distribution. We care about some summaries of the posterior

distribution.

Mm-hmm.

And then that informs when you use them. So I think that's the first

question. And

Mm-hmm.

I wanna emphasize that there's a lot of great work by Tamara Broderick

and her group, Justin Dumke and his colleagues about trying to get

more accurate variation inference. And sometimes that works really

well. Which is a bit of a... I can't really give you a more precise

prescription than that, because we have to go into the details of the

different problems.

But then the second point I made is sometimes you don't need a really

accurate answer. So what are examples of that? So in machine learning,

let's say you're just training a model, and maybe you're more interested

in either a more complicated model or using more data than improving the

accuracy of the inference. And then you

Mm-hmm.

look at something like performing a task, a classification prediction

and so forth under a computational budget,

Mm-hmm.

right? It turns out that it's better to have a sophisticated model with

very approximate inference than a less sophisticated model with more accurate

inference.

Mm-hmm.

Now, it's hard to know, and my big problem with variation inference

is it's hard to know which regime you're gonna be in. and to actually

justify it. And even once you've run variation inference, you don't have

that many diagnostics that can tell you, well, you know, you're doing

okay, but you would do much better if you improve the inference. So

I think that it's an open question. But then the other example that I

wanna bring up where we don't always need accurate inference is when

we're developing models. So this takes back into this idea of the Bayesian

workflow. Uh, that now has been championed. Uh, you know, so Andrew

Gellman and colleagues wrote a lot about it. Michael Bettencourt wrote

a lot about it. David Bly wrote about a lot about it. You know, arguably

George Box, right. Wrote a lot about it. And, um, and you know, if you,

if you ever work on an applied project, you sit down, you come up

with the first iteration of your model and, uh, arguably there are a

lot of problems with that model. And you don't need super accurate inference

to diagnose the problems with this model. You do a quick fit, a quick

approximation, and usually something obvious is gonna pop up.

Mm-hmm.

Then you revise the model. Then you do again, quick inference, okay? And

you keep refining and refining. And only once you have actually a

polished version of the model, do I think that it makes sense to, you

know, get out the big gun and the very accurate inference. And I

think that, you know, if we talk about pharmacometrics and epidemiology,

I'll give you some very precise examples of those situations.

Yeah, for sure. We'll do that in a few minutes. But yeah, thanks a lot for

that tour. That makes a lot of sense, actually, all of that. I will use

I realize my

those.

answers are very long, but your questions are very deep.

Yeah, yeah, yeah. No, no, that's really good. I mean, that's why the podcast

is for also, you know, going deep into the explanation that you cannot really

do in a paper, right? In the papers, usually it's a more technical audience

first.

Mm-hmm.

And so like, you're not going to really explain the difference between volitional

inference and MCMC in a paper, because if the audience already is supposed

to know that, well, why would you do that in a paper? So.

Yeah, and frankly the paper, so what ends up being the discussion

that ends up being in the paper is usually the discussion you've

had with the reviewer, which is a subset, which is really a subset

of everything you would like to discuss. So, yeah, it's nice to have

a more free format

Yeah.

to really think about those questions. And what I will say is that

actually the one format that was... that I really like where all these

questions come up is, you know, is teaching, is like workshops.

Yeah.

And so when you do a workshop on Bayesian modeling, you do a workshop on

Stan, on PMC or something like that, all the questions that I've brought

up, you know, how long should the sampling phase be? How long should

the warmup phase be? Should I use this or that algorithm? Those are questions

that a person sitting at a workshop Intro to Stan, intro to your

favorite language would ask. And so these end up being forums where

we do discuss these fundamental questions. Because even though they're deep,

they're elementary nonetheless. And I mean that in the most positive sense

of the word elementary possible.

Yeah, now completely you have completely outmasked the way I pick questions

for the show. I'm just doing the same as you did. Yeah, I teach a lot of

workshops too. And these questions are basically the questions that a lot of beginners

ask, where it's like they often have used variational inference because for

some reason, especially when they come from the classical machine learning

world, then using variational inference makes more sense to them because

it's closer from home, basically

Mm-hmm.

closer to home. And so yeah, like in the end the natural question is when

should I use variational inference? Why should I use it? Why should I even bother

with MCMC? Things like that. So I have a lot of the questions I've been

Yeah.

asking.

And these are totally open questions. I mean, correct me. Maybe you

have an answer that I missed. You know, we have heuristics and we

have good pointers and we have good case studies. But so much of this

remains unanswered. I think

Yeah, yeah,

not

yeah.

unapproachable. Let me be clear. Totally approachable. And you can,

you know, it's not crippling. We can totally still do a Bayesian modeling.

But there are open questions that linger.

Yeah. And then you have to strike the right balance between when you're answering

such a question from a beginner, where you want to be like intellectually honest

and saying,

Mm-hmm.

yeah, like these are still open questions, but at the same time, you don't

want them to walk away with the feeling that, well, this is completely...

like is completely undefined and I cannot really use these methods because

no, um, there is no clear rule of about what I would use in when

Mm-hmm.

and why. So, uh, like that's always a, an important balance to strike and

not, not always an easy one.

Yeah, yeah, I absolutely relate to that.

Yeah, so actually before we dig into a bit more of the applications, I'm

curious basically on the work you're doing right now because I really love

the fact that you're both working on MCMC and on approximate patient inference.

So I'm wondering what are the frontiers? currently in that field of algorithms

that you find particularly excited about. You already mentioned basically the

progress of the hardware, which opens a lot of avenues. I'm curious if there

are other things you have your eye on.

Yeah, I think that...

So. Well, I think hardware is an important question. And I think it's

a difficult question.

Mm-hmm.

And I wanna talk, I think, you know, maybe I wanna say a little bit more

because I think

Mm-hmm.

frontier is good. And not only that, I think it's an ambiguous frontier.

Mm-hmm.

Because I don't think, you know, to me, it's not clear. you know,

how much we're gonna get out of hardware for MCMC, for example,

and you know, what are gonna be the limits of that? And so what I'm

excited is that now we have GPUs, now we have several algorithms that

are GPU friendly. And I'll explain a little bit what that means,

but essentially, So at a very fundamental level, what it means is

you can run a lot of Markov chains

Mm-hmm.

in a synchronized time. And so you're not waiting for the slowest chain,

basically. That's kind of the intuition here. And people are like,

well, this is great. If you run a lot of chains, we can really make

the sampling phase much shorter. Like again, so let's say your target

effective sample size is 2000. And actually what effective sample

size you should target. That's another interesting and very fundamental

question. And it turns out a question where different people have

different opinions on

Mm-hmm.

within the Bayesian community. But let's say now for sake of argument, you're

targeting an effective sample size of 2000. And I'm not saying this

is what you should do or not do, but let's just say 2000. Then once

you run 2000 chains and you warm them up enough, right? But really

you only need one good sample per chain,

Thanks for watching!

which means your sampling phase can be one iteration. And then all your

questions about how long should the warmup be, it's no longer about

adapting the kernels so that I have a low auto correlation during the

sampling phase. Actually the auto correlation only matters in so far

as it reduces your bias. All right, and so

Mm-hmm.

I think that suddenly we've greatly simplified the question of diagnosing

how much computation we need to throw at the algorithm. If we have

a lot of variance, then suddenly it just becomes about bias decay.

And the sampling problem becomes much closer to the optimization problem.

Right. And then I can go into questions. So then there are interesting

things where you actually run a lot of chains. You can pool information

between the different chains to

Mm-hmm.

make the warmup phase shorter. And for some

Yeah.

problems, like you have these multimodal problems, like it's exponentially

shorter, like you're never going to get a good answer with a single chain

that doesn't use cross adaptation. Right. And so here, I want to give

a shout out to Marie-Lou Gabrielle and her work on MCMC that uses normalizing

flow for adaptation.

Mm-hmm.

Right. That's the thing, that's a technique where they actually, you

know, run 10,000 workers or chains. And I remember when I was talking

to her and colleagues, I was like, well, you know, once you've, you're

already running 10,000 chains to jump between the modes, you actually

don't need a long sampling phase at all,

Mm-hmm.

right? So that's one aspect of it. But then even for more ordinary problems,

you can show that the time it takes you to reduce the bias goes down

with the number of chains, because you are pooling information between

the Markov chains. And this is not something that we really understand.

Hmm.

Um, and so, you know, where I see the frontier is that actually,

if I run a lot of chains, I also get more accurate diagnostics. My computations

of our hat and its generalization become much more reliable. And I think the

Holy grail would be something where we don't have users specify, um, the

length of the warmup phase or the length of the sampling phase. We have

them think about what is your target ESS? That's the number of chains

that you run. And then we're going to automatically stop the warmup

phase. when we hit a certain target, right? And then suddenly,

we're starting to do optimal computation for MCMC. And I think that to do

optimal computation, at least in the way that I've described it, we

need those GPUs. And at the same time, I think that there are a lot

of problems that are not gonna be amiable to GPUs, right? It's still,

there's still this fundamental sequential component. which is the bias has

to go down, the warmup needs to happen, right? At some point, adding

more chains is not gonna help you. Whether the speed up you're gonna

get on this is not gonna be arbitrarily large, right? And then the benefit

you're gonna get from various reduction by running more chains, well,

once you've read. If your target ESS is 2000, maybe it doesn't help to

run 10,000 chains, right?

Yeah.

At least not immediately, right? So those are very clear, you know, questions

that arise about, you know, ultimately, what are, how far are we

going to be able to go with this ES, this, you know, running many

chains. But what I want to emphasize is that there's a, you know, there's

a computational game, which people think about, our algorithms are

faster. I also think there's a conceptual game. And there is the opportunity

to make MCMC more black box.

Mm-hmm.

We're less reliant on some very fundamental tuning algorithms.

And therefore, you know, we kind of get into this regime where the computation

that we're using is just the right amount. And that's really, you know,

if I have to say, like, I have another one or two years as a postdoc,

very optimistically, that's the problem I'd like to solve. Right? Like,

we remove the fundamental tuning parameters of MCMC. And there are other

approaches towards that. I'm not going to pretend that this is the only

angle to tackle this problem. Let me be absolutely clear. I think it's a

very, very promising.

Yeah, yeah,

Yeah,

for sure. That's fascinating.

what I will say is since we're going to go into applications is one

limitation

Mm-hmm.

with GPUs is they're very bad at solving ODEs.

Mm-hmm.

And a lot of the problems that I care about, you

Have

know,

a day!

have likelihoods. And in order to evaluate those likelihoods, you

need to solve an OD.

Mm-hmm.

So I don't think we're at a stage where running a thousand chains

on a GPU is going to solve all the problems in pharmacometrics that I

am, you know. deeply invested in, right? That said, you know, we have clusters

of CPUs where maybe we can run 60 to 120 chains and that can get us

some of the way.

Hmm, yeah. Yeah, I mean, when you were talking about that, I was thinking

that'd be super cool to have, you know, in Stan or Pimcey afterwards, like

at some point, optimize the way the number of chains and samples that are

taken instead of having them, because right now it's like, okay, we're gonna

run as many chains as we can with the GP or CPU we have. So that's

Mm-hmm.

already kind of automated. But the number of samples not really automated

is just a rule of thumb. Where it's like, OK, we think in general, these number

of samples work well. In PyMC, it's 1,000 per chain after warming up. But

what would be

Thanks

super

for watching!

cool is like, OK, PyMC or Stan, see what's there, the resources. And then it's

like, OK, so given the complexity of the posterior that we can see right now,

we are going

Yeah.

to run that many chains and that many sample per chain. That'd be super

cool, because also that would be something that would be easier for the beginners,

because they are

Mm-hmm.

really, really sometimes very anxious about having a lot of samples, even

though you know, no, you don't need 10,000 samples per chain for that simple

regression, but it's hard to explain.

Yeah, and so what ends up, also what ends up happening, very real

experience, and I do it myself, is the simple models that probably

don't need that many iterations. I run them for a lot of iterations because

that's computationally cheap to do. And then the hard models. where

each iteration is very expensive and the posterior distribution is

much more complicated, I actually end up, and where I would need more

iterations, I end up running less iterations.

Yeah.

Right, and I think a lot of people will sympathize with that. That's

my experience interacting with practitioners. I'll give you another example

of things I've seen my colleagues in epidemiology do is that when

their model starts getting really complicated, they start using somewhat

less over dispersed initializations.

Mm-hmm.

Even though we know that that's what we need for the convergence diagnostic

to be reliable. And I'll tell you a little bit more about that because

actually that's another question is like, what does it mean for initialization

to be over dispersed? And I have some answers to that are not quite,

you know, the regular answers. Um, but that's a huge problem. Like

when models get hard and again, those ODE based models. Right? You cannot

solve those ODE's if you throw insane parameter values at your model.

And so now people have to have to make compromises. And I think that,

you know, especially statisticians were a little bit cavalier. We try to, we tend

to be conservative, uh,

Mm-hmm.

because maybe in a way, you know, that's the role of the statistician

and the theorist is to be conservative and to be, and to play it safe.

But when the safe and the conservative heuristics, um become impractical

uh we need to think harder about you know okay if we don't want to be

too conservative but we still want to be safe what do we need to do

and that's where

Mm-hmm.

these questions of you know optimal warm-up length optimal sampling

phase optimal over dispersed initialization uh really come in into play

because if you're too conservative in your prescriptions You might

make your editor happy, but actually practitioners are going to have

a really hard time following those prescriptions.

Yeah,

And then they do

yeah.

things. I'm not saying that they do silly things instead, right? But

they do other things that are maybe less principled.

Mm-hmm. Yeah. Yeah, fascinating. I could spend the whole episode on these topics.

I really love it. Thanks a lot for diving so deep into these, Charles. But

let's get a bit more practical here and talk about

Uh huh.

what you do, basically, with epidemiology and pharmacometrics. So first,

can you define pharmacometrics for us? I personally don't know what that is,

Yeah.

how that differs from epidemiology and what do patient statistics bring to epidemiology

and pharmacometrics.

Yeah, so pharmacometrics, I mean, the way I would think about pharmacometrics

is pharmacometrics is to pharmacology, what econometrics is to economics.

Mm-hmm.

There's some people who want to emphasize that they're using quantitative

methods.

Mm-hmm.

Now, the particular field of pharmacometrics that I've worked on is called

PKPD modeling, so pharmacokinetics and pharmacodynamics. And essentially, what

happens is, let's say there's a new treatment that's being developed.

So it could be a new drug compound or it could be a number of things.

We're usually interested in two questions. One, how does the drug get

absorbed and how does it diffuse in the patient's body? That's called

the pharmacokinetics. And two, once the drug diffuses in the body,

what does it do to the body? And so that includes targeting a disease,

but also side effects. I was worried about how toxic the treatment

might be. And what people are trying to do with these models is based

on early data from clinical trials. So either on human beings or even on

animals, they try to predict. what is going to happen when we look at a broader

population of individuals and when we kind of like start changing the

treatment. So you have some, you know, some medical treatments can

get very complicated. You have questions of, you know, you know,

I have a certain drug compound. How much do I administer? How often

do I administer it? Is it better to take, you know, half a dose? every

half hour or only a single dose every hour. And so you have all these

combinatorics of possibilities. If I really increase the dose, do I immediately

get better effects? Or does it saturate? We often have these nonlinear

effects. These are called Michaelis-Menten models, where at some point,

adding more dose erases the cost and doesn't help the patient.

Yeah.

Right. And one way to do it would be

you run a ton of clinical trials and then you hope that something works

out that's extremely expensive and time consuming and, you know, well, Maybe

it's safer for the humans than the animals. Let me put it this way. Or

based on a bit of data, and then a really mechanistic model where you

actually really bake in some of your expertise as a pharmacologist,

as a biomedical engineer, you try to understand the underlying system

so that when you're trying out different regimens, you can really

Mm-hmm.

predict. what are gonna be the doses that are more promising? So it's

very useful to do exploratory analysis. And then it's also useful to do,

once a drug hits the market, you actually collect very imperfect data, you

actually have uncertainty, but you still want to keep learning about the

dose and the dosing regimen, right, from

Mm-hmm.

data from hospitals. And sometimes the data is rare, you have rare diseases,

rare conditions, and that kind of thing. And so that is... If you will,

then...

the domain of pharmacometrics. And what's extremely interesting is that within

that domain, you have, so first of all, you have what's called mechanistic

models. And what I mean by that is the parameters are interpretable.

The relationship between models

Mm-hmm.

is interpretable. Now, at contrast that with a neural network, for

example, where I might get good predictions, but then if I want to

do out of sample predictions, which is actually really what we want to do,

right? in pharmacometrics, right? Examples of other samples would be a

different dosing regimen, or I've tested the drug on an adult population.

What happens if I do it for children, right? That's the kind of pediatric

questions. We need to bake in the mechanistic understanding, which

doesn't exclude a role in neural networks. I think these can also

play a role, but I'll leave that for now. But then you have various

degrees details in the mechanism. You have some equations that are

very simple. So the two-compartment model with first-order absorption says the

human body is three compartments. There's the gut, where the drug arrives when

you already administer it. There is the central compartment where

the drug diffuses quickly. So usually that includes the blood. And then

maybe there's a peripheral compartment, so tissues where the drug diffuses

more slowly. That's obviously not a very detailed description of the human

body. And then you have models that are more complicated. So at Matron,

Matthew Riggs and some colleagues, they work on this bone mineral

density, where actually they had a lot of different parameters. And

now instead of having a system of three differential equations, you

have 30 differential equations. You have a ton of parameters, but you have

a lot of information, prior information, about these parameters.

Mm-hmm.

And then you have people who really throw differential equations with

you know, hundreds of states at them and, you know, thousands of

interpretable parameters. And frankly, I don't think we have the Bayesian

computation to fit those models, even though in theory, they lend themselves

extremely well to a Bayesian analysis, right? I think that realistically

we're somewhere in the semi-mechanistic regime. So these are models

that have some level of sophistication, but already we pay a dire price

for this sophistication, which is that the computation. can take hours

or days to fit. And so there's like really this potential for, better

Bayesian computation can really allow people to deploy better models

and more sophisticated. The other big aspect of pharmacometrics is

usually we have trials with data from different patients, there's

heterogeneity

Mm-hmm.

between the patients or similarity

Yeah.

between patients, so that lends itself very well to hierarchical modeling.

Yeah,

And

for

we

sure.

know hierarchical modeling is hard, right? It tends to create these posterior

distributions with geometries that are very frustrating. And I spent a lot

of time worrying about hierarchical models. I've done a lot of work

on the... nested Laplace approximations. So

Mm-hmm.

that's another nice example of an approximation. It's not variation inference

has very complimentary qualities. And what the nested Laplace approximation

allows you to do is marginalize out the latent variables in a hierarchical

model.

And people often explain nested Laplace as, oh, it's great, it reduces

the dimension of your problem. And then we can throw a quadratic

at the remaining parameters. We've had models where we had thousands of

hyperparameters, those were genetic problems where we were using horseshoe

parameters to select Gs. So even once you marginalize out, you know,

the latent variables, you still have a high dimensional problem.

So we threw Hamiltonian Monte Carlo at it, but it still made a big

difference because we simplified the geometry of the posterior distribution

by doing this marginalizations for the hierarchical models. So I'm very excited

about the prospect of having, you know, and that's the plus approximation

in Stan.

Thanks

We

for

have

watching!

a prototype. that works really well. We have some really cool automatic

differentiation supporting it.

Hmph.

But, and you know, the problem is again, I wanna try this on OD based

models.

Yeah.

I don't get a good approximation. I don't get efficient automatic differentiation.

I get something that's unstable. Now the simple examples where I got

it working, it actually gave surprisingly accurate results, right?

But this is again an example where here's this awesome. algorithm and statistical

methods, and it just gets frustrated by the nature of the problems we

encounter in pharmacometrics. Even though, you know, these are hierarchical

models and these methods are designed for hierarchical models. But again,

if your likelihood is not a general linear model, yeah, suddenly

those approximations become much more tricky. And that's why, that's why,

you know, I think that We have to use MCMC for these models.

Hmm.

Yeah, that's interesting. So in these cases, yeah, that makes it a bit clearer,

I think, for people in the practical cases where you would have to do that trade-off,

basically. How do you choose between the trade-off between the different

methods? I think it's very important.

And that said, I do want to say that there's some really cool approximations

that people

Mm-hmm.

do deploy in pharmacometrics.

Mm-hmm.

I recently read a...

I'm an anonymous reviewer, so I'm not going to give too many details.

Ha, sure.

But what I like, you know, it's, there were questions of what if we fix

those parameters or what if, you know, we draw these parameters from

their priors because they're removed away enough from the data,

Mm-hmm.

but maybe they're not influenced that much by the data, so the posterior

stays close to the prior, but nonetheless, we need that uncertainty for

the interpretable quantities in the middle. Right.

Hm-hm.

And so people are coming up with these compromises. And of course, now

we're again in the business of we have these awesome computational constraints.

People come up with these approximations either of the inference or even,

you know, they say, well, let's use a simpler model. We know there's

a more complicated model out there but maybe we still get all the answers

that we need with the simpler model. Right. And so now we get again

in this, you know, the question of understanding what are the simplifications

that we get away with? What are the ones where we pay a heavy price?

Can we actually quantify the price we're paying? Can we diagnose when

the simplification is too dire or not? And as far as I can tell,

the answer is, you know, we're not at the stage of diagnosing the problem,

but at least now people are taking this problem seriously enough that

they're building these case studies. Now based on these case studies

where, you know, we do try out the different models and we do fit

the complicated methods and we're able to say something about the simpler

methods because we did fit the complicated methods. So it's an academic

exercise in a way.

Mm-hmm.

But at least it gives us, you know, that's how it starts. And that's

how we're going to start developing the intuition and the heuristics

to then deal with, you know. deploying the approximations and the

simplifications. In practice.

Yeah, that makes sense. And so actually, do you have an example of a research

project that where you applied patient statistics in this field of pharmacometrics

and where really patient statistics helped you uncover very important insights?

thanks to all those methods that you've talked about.

Yeah. So I've never led a project in pharmacometrics. I've always collaborated

with pharmacologists.

Mm-hmm.

So, you know, and it's true that my work has been more methodological,

has been more developing towards in itself. That's it. What I can talk

about is some of the interactions I've had with pharmacometricians and some of

the works were maybe a more contributor than say, the project lead.

And I'll give you one example. And then I think that if we have time

and we talk a little bit about epidemiology, I have a very

Mm-hmm.

good example in epidemiology.

Ah, okay, I was...

But this is a, this was super cool actually. This is, this is a bit

anecdotal.

Mm-hmm.

I don't know that there's a preprint on this now, but I was in Paris,

I was visiting in CERN, which does a lot of medical work. And I was

visiting France Mautre's group, and they have fantastic people there.

And I'm talking with Julie Bertrand, and she's interested in pharmacogenetics,

right? So we're... We have some genome sequencing or some gene tracking

that's involved, and we're trying to see how a patient reacts to a treatment,

to a certain condition. And in that case, what ended up happening,

so one of the things that was frustrating is they're trying to identify

what are the genes that are meaningful, that seem to have a meaningful

connection to the outcome of the treatments. And they couldn't get anything

that is statistically significant in the traditional sense, right?

There wasn't one gene or one SNP that unembarrassed bigously stood

out as, yes, this is a meaningful snip and it should intervene in

oilier analysis.

Mm-hmm.

And what they had done was a Bayesian analysis where they had used

a horseshoe prior.

Mm-hmm.

And so what the horseshoe prior does is it does like, it's a regularization

tool that does a soft selection. And the way that soft selection

manifests is there is a quantity that you can look at. and it will

usually be bimodal. And one mode will indicate that the covariate that

corresponds to this snip is completely

regressed to zero. So that's an indication that this is not a very

useful explanatory variable. And then the second mode tells you actually

this variable matters and it's not regressed to zero. And so why

there are two modes is because there is uncertainty. But the very

cool thing is that even though there was no single snip that stood

out as this is the meaningful snip, uh, you had two snips that came up and

they both had a bimodal, right? So what that means is you couldn't definitely

say that, uh, either snip mattered, right? But you say, okay, with some

probability, either of the snips can matter.

Mm-hmm.

And here's where it gets very interesting is that actually, because

you have a multivariate posterior distribution, you can go a bit further

and you realize that the two snips are anti-correlated, right? And so

what that means is when you have a lot of posterior mass at one mode

for one snip that says, this variable matters, don't regress it to

zero, the other covariate would always get regressed to zero

Hmm.

and vice versa, right?

Mm-hmm.

So what the multivariate analysis tells you and what this proper treatment

of uncertainty tells you is like, yeah, you can't say that one snip

is statistically significant in the traditional sense, but now you have

this more comprehensive treatment of uncertainty that tells you, but

you know what? It has to be one of those two. You can't tell for

sure which one, but it has to be one of those two. And that's a nice

example where we're really, not just looking at the maximum likelihood

estimator or even the expectation value or just a variance, we're

really looking at essentially the quantiles, the extreme quantiles

Mm-hmm.

of the posterior distribution. And excuse me.

Sorry about that.

we're looking at multiple variables and

Mm-hmm.

the uncertainty across those multiple variables, right? So I think

that's a very neat example and you know. paper to look forward to when

it comes out or again. Again, this was anecdotal. This was a conversation

in the laboratory,

Yeah.

but I got very excited about that exam.

Yeah, I mean, that sounds super exciting. Thanks a lot for sharing that. And

yeah, like if when the paper is out, please get in touch and then that'd

be a fun thing to talk about again. And actually, so we're running a bit

long, but do you still have a few more minutes because I still have a few

questions for you. Or you want

Yeah,

to

yeah,

close out?

I can stick around for a bit. Yeah.

Okay, awesome. Yeah. So, uh, yeah, I'd like to, to talk to you about, uh,

priors and then maybe, um, Torsten, uh, or an epidemiology example, uh, or both

of them. Um, so yeah, maybe, um, basically you have, uh, been working on,

on Torsten, which is, if I understood correctly, uh, pharmacometrics application

of Stan models. So that was interesting to me because the Bayesian field

has been evolving really, really fast lately, especially with the new techniques,

the new software tools, and you've definitely been part of that effort,

Charles. So I'm wondering if any recent developments that have particularly

excited you, especially in your field of pharmacometrics and epidemiology,

and also in relationship with what you're doing with Tolston.

Yeah, okay, so let me talk about, let me give one example, right?

Mm-hmm, yeah.

So not a comprehensive answer, but an administrative answer.

Yeah, it's great.

Um, and so one, one fundamental, you know, yet another fundamental question

that would come up at a workshop, um, is okay. I have this ODE integrator,

right? I need to solve an ODE to evaluate my likelihood and the ODE's

come with certain tuning parameters. In

Mm-hmm.

particular, what is the precision with which you should solve your

ODE? And. And that is going to have an impact on the quality of your

entrance and also the run time. Right, because if you solve the OD

with a very strict tolerance, it takes much longer to solve that OD. Right.

And so again, that comes back to the question of how much computation

should we throw at the problem. And For the longest time, I didn't

have a good answer. And maybe I still don't to that question. And the

way this manifested is either, you know, at workshops when teaching

the subject, but even when we were writing the Stan manual, we have

a page on OD integrators. We are, we have a, you know, we state, these

are the default values that we use, but we don't have any clear recommendations

on what is the precision you should use. And we kind of assume that

the user knows their ODE well enough that they'll know which ODE integrator

to pick and what tolerance to set.

Mm-hmm.

which is not realistic. But it's just we didn't have an answer to that

question. And so recently, there was a paper. And so I know that.

Let me look up exactly who the... I know Aki is a co-author on it,

Aki the Tarii.

Mm-hmm.

But I want to give a shout out to the lead author.

Yeah, for sure. And we'll put that also in the show notes for this episode.

Yeah, so there you go. So Juho Timonen is the lead authors, is the

lead author, and then there are a bunch of other people on it. So

their paper is an important sampling approach for Bayesian ODIs. And so

essentially, what they realize is, if when you're solving, when you're

using a numerical integrator to evaluate your likelihood, really you're

not computing the true likelihood, you're computing an approximation

of this likelihood. And we have a lot of tools in statistics when we're

not dealing with the exact likelihood, but some approximation to this

likelihood, notably in the field of importance sampling.

Mm-hmm.

And what they came up with is a way to use those tools that exist

in importance sampling. to actually check whether the approximate likelihood

and therefore the tuning parameter of the ODE integrators are as precise

enough

Mm-hmm.

or not, right? And so that gives you a diagnostic that you can use.

It's not a completely perfect diagnostic and I think I'm still trying

to test that idea and play around with it and see how well it really

works. So I wanna try it out on pharmacometrics problems. Um, right

now I'm writing, uh, I've been tasked with writing, um, a tutorial on

Torsten. So we released part one a while ago. Now we're writing part

two and we promised that in part two, we'd explain how to tune the ODEs.

Mm-hmm.

Um, except we only know so much how to tune the ODEs. So I am trying

Yeah.

this method. I am trying this method, um, on the ODEs, but it's just

getting me, you know, thinking about. are the tolerances we're using

too conservative? Are they too strict? We actually get important

speed ups. I'm teaching a course this September in Leuven, the Advanced

Summer School in Basie Methods, where I'm gonna have the students

on an epidemiology problem, on a pharmacokinetic problem, try out different

tolerances and see the differences and then build this diagnostic

that's based on important sampling to check whether the precision with

which they're solving their ODEs is making meaningful changes to

the inference, right?

Mm-hmm.

And so again, I think that this is one tuning parameter where either

we're using ODE solvers that are not precise enough or ODE solvers

that are too slow. We're not being optimal in our computation. And this

is preventing us from either getting accurate answers or deploying

models with the sophistication that we would want to. And so that's a development

that I'm excited about. The one caveat that I will throw in this is

that right now we're still thinking about it as a single tuning parameter.

When what I've observed in practice is that the behavior of the ODE

can change widely depending on where you are in the parameter space.

Mm-hmm.

So for certain parameter values, you don't need to be super precise.

And for other parameter values, you need a lot more precision or you

need a different type of integrator because the OD just behaves in

a different way.

Yeah.

Now very concretely, how does this manifest? So I don't wanna say it's

hopeless, but what ends up happening is during the warmup phase,

where we start the Markov chains far off in the parameter space, or we

haven't tuned the MCMC sampler, so the Markov chains are still jumping

left and right.

Mm-hmm.

During the warm-up phase, we are more vulnerable to those extreme parameter

values. Whereas during the sampling phase, we get away with less strict

body solvers. So I think that somehow, what I would like to do is two

things. One, I would like to have a very automatic way of running this

diagnostics in Torsten. But I also wanna give users control over what

OD solver do they use at different stages of MCMC. They think that

makes a crucial difference. Another way to approach this problem is coming

up with good initializations. If I can start... my MCMC near, you know, within

the parameter space where I might land from the stationary distribution.

And I know that here, the parameter values are a bit less absurd.

And so solving the ODEs is a bit more feasible computationally. If

I can start there, then maybe I'm skipping the early regions that really

frustrate my ODE integrator and my MCMC sample. And the way this manifests

in practice is you'll have, you know, you run, let's say you run

eight chains. You have six of them that finish quickly, and then you

have two of them that are lagging because they're stuck somewhere

during the warmup phase. They're encountering this region where the

parameter values are a little bit absurd. Your OD is super hard to

solve, and that's eating up all your computation. And the truth

Mm-hmm.

is, you know, at least... And the way we do things right now, we always

wait for the slowest chain. By

Mm-hmm.

the way, we don't have to do it, right? So I'm excited about methods

to come up with good initializations. And I think that this is a place where variational

inference can be good, right? So the, especially

Aha.

now the pathfinder variational inference, right,

Mm-hmm.

was originally designed to produce good initializations for NCNC. And

so getting good initializations, that's a great example of, I need a good answer,

but not a super precise answer. Right. Uh, and so, you know, if somehow

pathfinder can help me skip the regions that frustrate OD integrators,

I think that's a big win for pharmacometrics. Again, that's, that's

something we have to test and really try out. But I will say now going

back, I'm going to make another connection back to our hat. which is

when we think about over dispersion. really what over dispersion means.

So we have shown a bit formally that what makes our hat reliable and

we define reliability in a formal sense is the initial variance has

to be large relative to the initial bias. Now, if you have an initialization,

like you draw your sample from your prior and that's reliable,

Mm-hmm.

and then you throw variation inference, right? If variation inference reduces

your squared bias more than it reduces your variance, it turns out

you preserve the property of reliability. So there's actually a sense

that we might be able to get good initializations for MCMC without

compromising, uh, the reliability of our diagnostics for convergence,

right? And these are all the pieces that I think can come together and

really help a great deal with, you know, pharmacometrics, but more generally,

uh, with OD based models. and even more generally, models based on implicit

functions. And I do include things that use nested Laplace approximations

in that, because that's an optimization problem, that's an implicit

function. It has the same kind of misbehaviors that an ODE has, but also

technologies that we develop for ODE is in a Bayesian context, better initializations,

different tolerances, important sampling corrections would apply

to nested Laplace stuff. So those are the things that I'm excited about

Um, but it's going to take time.

Hey,

I just

yeah.

want to, I just want to be perfectly honest. It takes time to really,

you know, uh, you know, between the paper and. The software implementation

and, uh, you know, clear description in the manual that the users can,

can follow. It takes a lot of time.

Yeah, yeah, for sure. I mean, that stuff is really at the frontier of the

research. So it does make sense that it takes time to percolate

Yeah, yeah.

from how to find a solution to, okay, this is how we can implement it and

reproduce it reliably.

Yeah, exactly. And I'm most concerned about how speculative I am for

some of the ideas I'm sharing. But I think these are directions where

it's worth pushing the research. That has the potential to have a

really big impact.

Yeah, yeah, for sure, for sure. And so I put in the show notes the Pathfinder

paper, actually. That made me think that I should do an episode about the

Pathfinder paper, basically, what that is about and what that means concretely.

Mm-hmm.

So yeah, I'll try to do that

Yeah.

in

And

the

if you haven't

near

reached

future.

out to Lu Zhang, you know, I mean, it, you know, first of all, the paper

is great. But actually sitting down and discussing this with her,

you know, at the blackboard or whatever, like, she has so many ideas that

have not appeared in the paper itself. I think, you know, if I can recommend

a guest. If you haven't had her already, I don't know if you... But

yeah, Lu

No,

Zheng, I think,

I didn't.

would be fantastic to interview.

Yeah, yeah, exactly. I was actually thinking about inviting her on the podcast

to talk about Pathfinder because she's the lead author on the paper and also

all the other authors have been on the show, Bob

Ah

Carpenter,

cool.

Andrew

Cool.

Gellman, Akive Tareiso. Lu Zheng is missing, so definitely

Mm-hmm.

need to correct that. So yeah, in the near future, definitely try. to have

that episode that would be a very interesting one. So maybe can you talk

Charles about... So I'll give you two avenues and you pick the one you prefer

because I don't want to take too much of your time but basically I'm curious

either about hearing an example of your work in the epidemiology field Or you

can talk a bit more in general about a very, very common question that students

always ask me, and it's about priors. So basically, how do you choose the

prior, how do you approach the challenge of choosing appropriate prior distributions,

and especially when you're dealing with complex models. So these are the two avenues

I have in mind. And feel

Okay.

free to... Pick one or... Pick both.

Okay, so let me say that about priors, I think that...

I want to get better answers,

Yeah.

a better answer to this question. And, you know, like the next time

workshop that I'm giving, I don't have a module on priors that I find

satisfactory. And so I'm still undergoing this journey. But what I will

do is I'll give, I'll talk about epidemiology and I will talk about

the prior that we use there. So that will be, that would be an example.

And, you know, I like to think. through examples, I like to think

Mm-hmm.

through the anecdotal, as complementary to the formal, right? I'm

a big fan of fairy tales and fables, simple stories that have good themes.

Yeah.

So what happened in epidemiology, and this will be a good example, is, well,

the pandemic happened, and suddenly COVID, we're all going home.

And we had colleagues in epidemiology, in particular Julien Rioux,

who actually met Julien Rioux at the StankCon. He was in Cambridge,

he was a PhD student at the time and he demonstrated his model and

he was using those ODI based models. So now instead of having a drug

compound that flows between different parts of the body, what you have

is... you separate the population, the human population, to these

different compartments. So susceptible individuals, infected individuals,

recovered individuals, and then the individuals flow between the compartments.

And there are a bit more layers. But basically the mathematical formalism

that I had familiarized myself with in the context of pharmacometrics

turned out to be very relevant to... certain classes of epidemiological

models. And essentially, Julien was working on an early model of COVID-19.

They were trying to estimate the mortality rate. There was a lot of uncertainty

in the data. There were a lot of

things to correct for. Right. So, for example, early on, not everyone

got tested and testing was not widely available.

Mm-hmm.

Who got tested? The people with severe symptoms. So now if you think about

you're trying to estimate the mortality rate

Yeah.

and according to your data, the people who catch the disease are

only the ones who have severe symptoms, then it looks like a lot of

people are dying from the disease. I mean, a lot of people were dying

from the disease, but that inflates the number. There's a bias

Yeah,

because you're

yeah,

only

for sure.

testing the people who are sick. You're not testing the people with

mild symptoms or even with no symptoms.

No, for sure. So, clear sampling by ASEA.

The other bias was some of the people who had caught the virus had

not died yet. So you count them as living, but that doesn't mean they

survived the disease. So that's a bias

Yeah.

in another direction. And that's an example where you actually have

a somewhat mechanistic model. That's based on the epidemiology of

how a disease transmits and circulates in a population. Then on top

of that, you need to build a measurement model to account for, you know,

how the data is collected. Right. But at the end of the day, none of the,

you know, we were not able to draw any conclusions. unless we understood

what was the rate of people who were symptomatic or had severe symptoms.

Right. And so there's one parameter in the model, which is the asymptomatic

rate. And so now you have two options in a classic classical statistics

framework. Either you fix the parameter and then you're making a

strong assumption and maybe you try different values of the fixed

parameter. Right. Or you just say, well, I don't know this. So really

I have no idea what the mortality rate is because maybe the entire

population was infected or maybe only a small fraction was infected

and everyone in that small fraction had the severe disease, right? And

so we needed an in-between, between saying we don't know anything and saying

we know everything. And this is why I think we're basing shines, which

is we can quantify uncertainty. We have more nuanced statements. through the

language of probability about what our state of knowledge is. And actually

what had happened is there were some instances where we had measured asymptomatic

rates. The example was a cruise ship, the Diamond Princess, in the

coast of Japan. So they had identified some cases of COVID-19. They

put the cruise ship in quarantine and they test everybody, regardless

of whether they had symptoms or not. So now you have a small population

and based on that small population, you get an estimate of the

Yeah.

asymptomatic rate. And then you had one or two incidents where people

do that. There were some cities where they had done some other experiments

and measured some other data. And so then you bring all that information

and you use that information. to construct priors and then to propagate uncertainty

into the model, right? The reason we're able to make predictions and

then to

calculate things like the mortality rate with the appropriate uncertainty

is because we had a prior on the asymptomatic rate. And that's a very

nice example. This is more an example of why it was crucial to have

a prior rather than, you know. How should you construct priors in general,

right? This is a bit of a specific case,

Yep.

but it's a very good example. And so I'll recommend two papers. One is the

one by Julien Rioux, which was, so this was disappeared in PLOS Medicine.

And it's with a lot of contributors, but it's an estimation of seroscoped

tumor mortality during the early stages of an epidemic. And then the

other paper, that goes a bit more into what are the lessons that we

learned from this, from a Bayesian workflow perspective, is Bayesian

workflow for disease transmission modeling in Stan. And so here the first author

was Leo Grintjerdsch.

Mm-hmm.

And then we had also Lisa Semenova and Julien Rioux as co-authors.

Mm-hmm.

And this is a beautiful paper that really goes into the, you know,

here's the prior reuse. Here's the first version of the model. Here

are the limitations with this model. Here's how we diagnose the models.

Here's the next iterations. And we go through all the iterations. And I

like to show the numbers that the model that eventually we used to model

COVID-19 was the 15th iteration. And along the way we had the model

that took three days to run. And we had to change the way we wrote

the Stan model to improve the computation. So that had to do with how

the OD was parameterized and how the automatic differentiation was happening.

We got it from three days to two days. And that's useful not, sorry,

three days to two hours, not two days, two hours, drastic speed up.

And so that not only was good because the inference is faster, but that's

what allowed us to then use more sophisticated versions of the model.

And so

Yeah.

all that is described in that. Bayesian workflow for disease transmission.

paper.

Hmm, oh yeah, need

So

that.

that's my epidemiology fairy tale. And I don't mean

Yeah.

that in the sense that it had a happy ending or that everything was

nice and glowing. I mean that this is a very nice story that touches

upon a lot of interesting things.

Yeah, yeah, for sure. Definitely going to use link to this paper in the show

notes. Already have it on the Stan website. Actually, you've done a case

study on this, so that's perfect with all the co-authors. So I'm going to put

that right now in the show notes. And maybe last question before letting you

go, Charles. I'm breaking records these days on the episodes. Like episode 89

is going out this week actually. And it's so far the longest episode. Uh, it's

about two hours. And right now we're like approaching this record, uh, Charles.

So, uh, like well done. And at the same time.

Okay. There's a fantastic podcast and I forget the name. It's a German

podcast.

Uh huh.

I can't believe I forget the name. I'll try and email you what it is,

but basically they do interviews with these, you know, intellectuals

and well established, you know, politicians and well, I'm not saying

all politicians are intellectuals, but people have opinions and thoughts and there's

no time limit to the interview. And they just go on and on and on and on and

on. And they just have so many topics to discuss. It's really, I'm

not saying you should do that with me, but you could imagine like,

you know, someone like, you know, Some of the other co-authors that

have come up, I feel like you could talk six hours with them and you

would just

Oh

pick

yeah,

their

for

brand

sure.

and they would have so much insights to

Yeah.

share.

No, for sure.

So

Like, I mean, most

they

of

are.

the time the limitation is the, is the guest's own, own time.

Right, right,

Yeah.

right. But the

Yeah, for

dev,

sure.

yeah, I think it's really cool, this idea of, yeah, we're gonna do, what

if we didn't have a time limit on the interview? What would happen?

And

Exactly.

eventually somebody gets hungry and that's what happens, but. Ha ha

ha.

Yeah, yeah, exactly. Yeah, but so basically before asking you the last two

questions, I'm also wondering, because you also teach, so that's super interesting

to me. And I often hear that a lot of practitioners and beginners in particular

might be hesitant or even intimidated. to adopt patient methods because they perceive

them as complex. So I don't necessarily agree with that underlying assumption,

but without necessarily disputing it, what do you do in those cases? Basically,

what resources or strategies do you recommend to those who want to learn

and apply patient techniques in the work, but might be intimidated or hesitant?

Yeah, okay, so I think that usually when I teach.

Most of the time it's people who are already interested in Bayesian

methods,

Mm-hmm.

especially if it's a workshop on stand, the people who do sign up, they already

have adhered to Bayesian methodology. So I don't know that I've really

had to convince that many people. Definitely when I was a TA at Columbia,

Mm-hmm.

I had those conversations a bit more with the

Mm-hmm.

students, especially when I TA'd the PhD level. that everyone has to

take applied statistics and

Hm.

not everyone is gonna do Bayesian, right? And so we have those conversations.

But I think that ultimately what it is, is not that Bayesian is complex,

it's analysis is complex. Analysis is difficult. And when you use less

complicated methods,

like, maximum likelihood estimates or point estimates. And these are

not simple. I really don't want to undermine the difficulty related to

those methods and the fact that they are useful in a lot of applications.

But anyway, those methods are simple and they will work for simple analysis.

But if an analysis requires you quantify uncertainty, to propagate uncertainty,

to take some decisions with imperfect data. And then someone says,

well, I don't wanna be Bayesian because that's complicated, but I

still want to quantify uncertainty and I still want to propagate uncertainty

and I still want to make predictions. Well, suddenly the classical methods,

you really have to do a lot of gymnastic to get them to do what you're

interested in. And And so the classical methods also become complicated.

And that's more because you're trying to use them for, uh, you know,

to do a sophisticated analysis. And so I think that what matters

is, you know, to really meet practitioners, uh, where they are, what

is the problem they're interested in? And. you know, does the, you

know, is it, is it a problem where, you know, the, the complexity of the

analysis would be handled in a relatively straightforward way by a Bayesian

analysis? And the answer could be yes, it could also be no. And if it's

no, that's fine. You know, I think that again, we have, we have to be

loyal to, to the problem, more so than to the field or to, or to the

method, but that's how I would start this conversation. And so when,

you know, when I'm sitting with pharmacometricians, And by the way, a

lot of them don't use Bayesian. I think I can talk for 30 minutes

about what do we want to get out of a pharmacokinetic analysis without

bringing up what kind of statistics we do. And once we've established

the common goals, then we think about what are the methods that are

going to get there. That's how I would do it. Frankly, I haven't had that

much experience doing it. So take what I say with a grain of salt.

Yeah, but I really think it's, Bayesian is complicated because it confronts

you with the complexity of data analysis. Right? It doesn't

Yeah.

introduce the complexity. Let me put it this way.

Yeah, yeah. No, completely agree. And that usually is... I usually give an answer

along those lines. So

Mm-hmm.

that's interesting. To see that we have converged, even though we didn't, we

didn't, you know, talk about that before the show. so that people don't

start having a conspiracy theories about Charles and I trying to push back the

message.

This is not rehearsed.

Exactly. Awesome. Well, is there a topic I didn't ask you about that you'd like

to mention before we close up the show?

Um, no. Let me keep

Perfect.

it simple.

Yeah. I mean, we've talked about a lot

There

of

are

things.

a lot of topics that I think we could have very interesting conversations

about, but I also think that we're starting to hit the time constraints.

Awesome. Well, let's go to show then. Me too, I could keep asking you a lot

of things, but let's do another episode. Another day. That'd be a fun thing.

Or if one day you have a model you've been working about and that you think

would be beneficial for listeners to go through. I have this new format now,

which are the modeling webinars, where basically you would come on the show

and share your screen and show us your code and do some live coding about

basically a project that you've been working on. And so if one day you have

something like that, feel free to get in touch and we'll get that organized

because... It's a really cool new format and I really love it. In August,

we've had Justin Boyce showcasing how to do Bayesian modeling workflow in the

biostatistics world. So even if you're not a biostatistician, it's a really

useful thing because basically, as we were saying, Bayes are basically methods.

Even though the field might... not be yours, the methods are definitely transferable

in the workflow. So that was a very fun one. And then we've got one coming

up in September in theory with Benjamin Vincent and we're going to dive into

the new do operator that we have in Pimc and how to use that. So that's

going to be a very fun one too.

Yeah, I love the idea of that format, and I'm definitely going to

check out the two webinars.

Yeah,

The

yeah.

one that's already there and the one that's coming up. Actually, the

two topics sound extremely interesting.

Yeah, I'll send that to you. And yeah, for sure that that's something I've

been starting to do. And so if you have one such analysis one day, feel

free to reach out that'd be a fun one. Um,

Absolutely.

so before letting you go, of course, I'm going to ask you the last two questions

I ask every guest at the end of the show. One, if you had unlimited time

and resources, which problem would you try to solve?

Oh yeah, I thought about this and

Hehe.

I feel like so much of my work is about working under computational

constraints,

Mm-hmm.

you know, and time constraints and resource constraints and suddenly you relax

all of this. And...

I think. For me, they're really big questions.

let's say at a curiosity level, right? I think there are, you know,

there are utilitarian questions and philanthropic questions that

I might prioritize. But my thinking right now is, I really love the

problems I worked on in astronomy. I am mind blown by some of the stuff

my colleagues do in cosmology, where they're trying to understand, you know,

the early universe. They have models with six parameters that apparently