Reproducible
Statistical Inference
Jonathan Huggins
Department of Mathematics & Statistics+ Faculty of Computing & Data Sciences
Boston University
Based on work with Jeﬀ Miller and Jiawei Li

1

Bayesian inference2Bayesian inference
interest 𝜃 [e.g., “skill” of a baseball player]• Goal: learn about unobserved phenomenon (parameter) of2

Bayesian inference
interest 𝜃 [e.g., “skill” of a baseball player]• Goal: learn about unobserved phenomenon (parameter) of• Prior beliefs π0(θ) about the phenomenon

prior

✓
<latexit sha1_base64="sKlV6FZOXFp9/WXwNNdnekD6QqY=">AAACBHicbVDLTgJBEJzFF+IL9ehlIjHxItlFEz0SvXjERB4JbMjsMLADM7ObmV4TsuHq2at+gzfj1f/wE/wLB9iDgJV0UqnqTndXEAtuwHW/ndza+sbmVn67sLO7t39QPDxqmCjRlNVpJCLdCohhgitWBw6CtWLNiAwEawaju6nffGLa8Eg9wjhmviQDxfucErBSowMhA9ItltyyOwNeJV5GSihDrVv86fQimkimgApiTNtzY/BTooFTwSaFTmJYTOiIDFjbUkUkM346u3aCz6zSw/1I21KAZ+rfiZRIY8YysJ2SQGiWvan4n9dOoH/jp1zFCTBF54v6icAQ4enruMc1oyDGlhCqub0V05BoQsEGtLBlGIYXkhtamRRsNt5yEqukUSl7l+XKw1WpepullEcn6BSdIw9doyq6RzVURxQN0Qt6RW/Os/PufDif89ack80cowU4X7/hmJhy</latexit>

2

Bayesian inference
interest 𝜃 [e.g., “skill” of a baseball player]• Goal: learn about unobserved phenomenon (parameter) of• Prior beliefs π0(θ) about the phenomenon•

Observe data x [e.g., player’s performance in somegames]

prior
likelihood✓<latexit sha1_base64="sKlV6FZOXFp9/WXwNNdnekD6QqY=">AAACBHicbVDLTgJBEJzFF+IL9ehlIjHxItlFEz0SvXjERB4JbMjsMLADM7ObmV4TsuHq2at+gzfj1f/wE/wLB9iDgJV0UqnqTndXEAtuwHW/ndza+sbmVn67sLO7t39QPDxqmCjRlNVpJCLdCohhgitWBw6CtWLNiAwEawaju6nffGLa8Eg9wjhmviQDxfucErBSowMhA9ItltyyOwNeJV5GSihDrVv86fQimkimgApiTNtzY/BTooFTwSaFTmJYTOiIDFjbUkUkM346u3aCz6zSw/1I21KAZ+rfiZRIY8YysJ2SQGiWvan4n9dOoH/jp1zFCTBF54v6icAQ4enruMc1oyDGlhCqub0V05BoQsEGtLBlGIYXkhtamRRsNt5yEqukUSl7l+XKw1WpepullEcn6BSdIw9doyq6RzVURxQN0Qt6RW/Os/PufDif89ack80cowU4X7/hmJhy</latexit>

2

Bayesian inference
interest 𝜃 [e.g., “skill” of a baseball player]• Goal: learn about unobserved phenomenon (parameter) of• Prior beliefs π0(θ) about the phenomenon•

Observe data x [e.g., player’s performance in somegames]•

Assume a probabilistic model p(x ∣ θ)

priorposterior
likelihood✓<latexit sha1_base64="sKlV6FZOXFp9/WXwNNdnekD6QqY=">AAACBHicbVDLTgJBEJzFF+IL9ehlIjHxItlFEz0SvXjERB4JbMjsMLADM7ObmV4TsuHq2at+gzfj1f/wE/wLB9iDgJV0UqnqTndXEAtuwHW/ndza+sbmVn67sLO7t39QPDxqmCjRlNVpJCLdCohhgitWBw6CtWLNiAwEawaju6nffGLa8Eg9wjhmviQDxfucErBSowMhA9ItltyyOwNeJV5GSihDrVv86fQimkimgApiTNtzY/BTooFTwSaFTmJYTOiIDFjbUkUkM346u3aCz6zSw/1I21KAZ+rfiZRIY8YysJ2SQGiWvan4n9dOoH/jp1zFCTBF54v6icAQ4enruMc1oyDGlhCqub0V05BoQsEGtLBlGIYXkhtamRRsNt5yEqukUSl7l+XKw1WpepullEcn6BSdIw9doyq6RzVURxQN0Qt6RW/Os/PufDif89ack80cowU4X7/hmJhy</latexit>

2

Bayesian inference
interest 𝜃 [e.g., “skill” of a baseball player]• Goal: learn about unobserved phenomenon (parameter) of• Prior beliefs π0(θ) about the phenomenon•

Observe data x [e.g., player’s performance in somegames]

prior
posterior
likelihood

Assume a probabilistic model p(x ∣ θ)•
• Combine prior & likelihood to form posterior:p(x ∣ θ)π0(θ)π(θ ∣ x) =p(x)✓<latexit sha1_base64="sKlV6FZOXFp9/WXwNNdnekD6QqY=">AAACBHicbVDLTgJBEJzFF+IL9ehlIjHxItlFEz0SvXjERB4JbMjsMLADM7ObmV4TsuHq2at+gzfj1f/wE/wLB9iDgJV0UqnqTndXEAtuwHW/ndza+sbmVn67sLO7t39QPDxqmCjRlNVpJCLdCohhgitWBw6CtWLNiAwEawaju6nffGLa8Eg9wjhmviQDxfucErBSowMhA9ItltyyOwNeJV5GSihDrVv86fQimkimgApiTNtzY/BTooFTwSaFTmJYTOiIDFjbUkUkM346u3aCz6zSw/1I21KAZ+rfiZRIY8YysJ2SQGiWvan4n9dOoH/jp1zFCTBF54v6icAQ4enruMc1oyDGlhCqub0V05BoQsEGtLBlGIYXkhtamRRsNt5yEqukUSl7l+XKw1WpepullEcn6BSdIw9doyq6RzVURxQN0Qt6RW/Os/PufDif89ack80cowU4X7/hmJhy</latexit>

2

Bayesian inference
interest 𝜃 [e.g., “skill” of a baseball player]• Goal: learn about unobserved phenomenon (parameter) of• Prior beliefs π0(θ) about the phenomenon•

Observe data x [e.g., player’s performance in somegames]

prior
posterior
likelihood

Assume a probabilistic model p(x ∣ θ)•
• Combine prior & likelihood to form posterior:p(x ∣ θ)π0(θ)π(θ ∣ x) =p(x)✓<latexit sha1_base64="sKlV6FZOXFp9/WXwNNdnekD6QqY=">AAACBHicbVDLTgJBEJzFF+IL9ehlIjHxItlFEz0SvXjERB4JbMjsMLADM7ObmV4TsuHq2at+gzfj1f/wE/wLB9iDgJV0UqnqTndXEAtuwHW/ndza+sbmVn67sLO7t39QPDxqmCjRlNVpJCLdCohhgitWBw6CtWLNiAwEawaju6nffGLa8Eg9wjhmviQDxfucErBSowMhA9ItltyyOwNeJV5GSihDrVv86fQimkimgApiTNtzY/BTooFTwSaFTmJYTOiIDFjbUkUkM346u3aCz6zSw/1I21KAZ+rfiZRIY8YysJ2SQGiWvan4n9dOoH/jp1zFCTBF54v6icAQ4enruMc1oyDGlhCqub0V05BoQsEGtLBlGIYXkhtamRRsNt5yEqukUSl7l+XKw1WpepullEcn6BSdIw9doyq6RzVURxQN0Qt6RW/Os/PufDif89ack80cowU4X7/hmJhy</latexit>• Benefits: coherent belief updates, uncertainty
quantification, flexible modeling, and more2

Bayesian inference
interest 𝜃 [e.g., “skill” of a baseball player]• Goal: learn about unobserved phenomenon (parameter) of• Prior beliefs π0(θ) about the phenomenon•

Observe data x [e.g., player’s performance in somegames]

prior
posterior
likelihood

Assume a probabilistic model p(x ∣ θ)•
• Combine prior & likelihood to form posterior:p(x ∣ θ)π0(θ)π(θ ∣ x) =p(x)✓<latexit sha1_base64="sKlV6FZOXFp9/WXwNNdnekD6QqY=">AAACBHicbVDLTgJBEJzFF+IL9ehlIjHxItlFEz0SvXjERB4JbMjsMLADM7ObmV4TsuHq2at+gzfj1f/wE/wLB9iDgJV0UqnqTndXEAtuwHW/ndza+sbmVn67sLO7t39QPDxqmCjRlNVpJCLdCohhgitWBw6CtWLNiAwEawaju6nffGLa8Eg9wjhmviQDxfucErBSowMhA9ItltyyOwNeJV5GSihDrVv86fQimkimgApiTNtzY/BTooFTwSaFTmJYTOiIDFjbUkUkM346u3aCz6zSw/1I21KAZ+rfiZRIY8YysJ2SQGiWvan4n9dOoH/jp1zFCTBF54v6icAQ4enruMc1oyDGlhCqub0V05BoQsEGtLBlGIYXkhtamRRsNt5yEqukUSl7l+XKw1WpepullEcn6BSdIw9doyq6RzVURxQN0Qt6RW/Os/PufDif89ack80cowU4X7/hmJhy</latexit>• Benefits: coherent belief updates, uncertainty
quantification, flexible modeling, and more

• Assumption: measurement model correct: observed x has
distribution p(x ∣ θtrue)2Motivating example:
Bayesian model selection3Motivating example:
Bayesian model selection• Given data x, select betweena (finite or countable) set of models

m1, m2, …

3

Motivating example:
Bayesian model selection• Given data x, select betweena (finite or countable) set of models

m1, m2, …

• Running example: systematics3Motivating example:
Bayesian model selection• Given data x, select betweena (finite or countable) set of models

m1, m2, …

m1

Fin whale<latexit sha1_base64="psq2jvV0jDNdqjwPulyjlu1ZWqw=">AAACAXicbVDLSgNBEOz1GeMr6tHLYBC8GHajoMeAF48RzQOSJcxOZrNjZmaXmVkhLDl59qrf4E28+iV+gn/hJNmDSSxoKKq66e4KEs60cd1vZ2V1bX1js7BV3N7Z3dsvHRw2dZwqQhsk5rFqB1hTziRtGGY4bSeKYhFw2gqGNxO/9USVZrF8MKOE+gIPJAsZwcZK96Ln9Uplt+JOgZaJl5My5Kj3Sj/dfkxSQaUhHGvd8dzE+BlWhhFOx8VuqmmCyRAPaMdSiQXVfjY9dYxOrdJHYaxsSYOm6t+JDAutRyKwnQKbSC96E/E/r5Oa8NrPmExSQyWZLQpTjkyMJn+jPlOUGD6yBBPF7K2IRFhhYmw6c1seo+hcME2q46LNxltMYpk0qxXvolK9uyzXanlKBTiGEzgDD66gBrdQhwYQGMALvMKb8+y8Ox/O56x1xclnjmAOztcvFi2W3Q==</latexit>

Blue whale
Grey whale• Running example: systematics• Goal: learn about evolutionary history of a

Minke whale

group of species [e.g., whales]m2

Fin whale<latexit sha1_base64="fhdmpccEdfAZWx7deDX0E4BpTnA=">AAACAXicbVDLSgNBEOz1GeMr6tHLYBC8GHajoMeAF48RzQOSJcxOZrNjZmaXmVkhLDl59qrf4E28+iV+gn/hJNmDSSxoKKq66e4KEs60cd1vZ2V1bX1js7BV3N7Z3dsvHRw2dZwqQhsk5rFqB1hTziRtGGY4bSeKYhFw2gqGNxO/9USVZrF8MKOE+gIPJAsZwcZK96JX7ZXKbsWdAi0TLydlyFHvlX66/ZikgkpDONa647mJ8TOsDCOcjovdVNMEkyEe0I6lEguq/Wx66hidWqWPwljZkgZN1b8TGRZaj0RgOwU2kV70JuJ/Xic14bWfMZmkhkoyWxSmHJkYTf5GfaYoMXxkCSaK2VsRibDCxNh05rY8RtG5YJpUx0WbjbeYxDJpViveRaV6d1mu1fKUCnAMJ3AGHlxBDW6hDg0gMIAXeIU359l5dz6cz1nripPPHMEcnK9fF8iW3g==</latexit>

Blue whale
Grey whale
Minke whale

.
.
.
3

Motivating example:
Bayesian model selection• Given data x, select betweena (finite or countable) set of models

m1, m2, …

m1

Blue whale
Grey whale
Minke whale

group of species [e.g., whales]consistent with observed species
characteristics x [e.g., genetic data, physicalfeatures such as coloring and size]p(x ∣ mi)π0(mi)π(mi ∣ x) =∑j p(x ∣ mj)π0(mj)Fin whale<latexit sha1_base64="psq2jvV0jDNdqjwPulyjlu1ZWqw=">AAACAXicbVDLSgNBEOz1GeMr6tHLYBC8GHajoMeAF48RzQOSJcxOZrNjZmaXmVkhLDl59qrf4E28+iV+gn/hJNmDSSxoKKq66e4KEs60cd1vZ2V1bX1js7BV3N7Z3dsvHRw2dZwqQhsk5rFqB1hTziRtGGY4bSeKYhFw2gqGNxO/9USVZrF8MKOE+gIPJAsZwcZK96Ln9Uplt+JOgZaJl5My5Kj3Sj/dfkxSQaUhHGvd8dzE+BlWhhFOx8VuqmmCyRAPaMdSiQXVfjY9dYxOrdJHYaxsSYOm6t+JDAutRyKwnQKbSC96E/E/r5Oa8NrPmExSQyWZLQpTjkyMJn+jPlOUGD6yBBPF7K2IRFhhYmw6c1seo+hcME2q46LNxltMYpk0qxXvolK9uyzXanlKBTiGEzgDD66gBrdQhwYQGMALvMKb8+y8Ox/O56x1xclnjmAOztcvFi2W3Q==</latexit>• Running example: systematics• Goal: learn about evolutionary history of a
• Approach: infer which phylogenetic trees areπ(m1 ∣ x) = 0.8m2

π(m2 ∣ x) = 0.1

Fin whale<latexit sha1_base64="fhdmpccEdfAZWx7deDX0E4BpTnA=">AAACAXicbVDLSgNBEOz1GeMr6tHLYBC8GHajoMeAF48RzQOSJcxOZrNjZmaXmVkhLDl59qrf4E28+iV+gn/hJNmDSSxoKKq66e4KEs60cd1vZ2V1bX1js7BV3N7Z3dsvHRw2dZwqQhsk5rFqB1hTziRtGGY4bSeKYhFw2gqGNxO/9USVZrF8MKOE+gIPJAsZwcZK96JX7ZXKbsWdAi0TLydlyFHvlX66/ZikgkpDONa647mJ8TOsDCOcjovdVNMEkyEe0I6lEguq/Wx66hidWqWPwljZkgZN1b8TGRZaj0RgOwU2kV70JuJ/Xic14bWfMZmkhkoyWxSmHJkYTf5GfaYoMXxkCSaK2VsRibDCxNh05rY8RtG5YJpUx0WbjbeYxDJpViveRaV6d1mu1fKUCnAMJ3AGHlxBDW6hDg0gMIAXeIU359l5dz6cz1nripPPHMEcnK9fF8iW3g==</latexit>

Blue whale
Grey whale
Minke whale

.
.
.
3

Motivating example:
Bayesian model selection• Given data x, select betweena (finite or countable) set of models

m1, m2, …

m1

Blue whale
Grey whale
Minke whale

group of species [e.g., whales]consistent with observed species
characteristics x [e.g., genetic data, physicalfeatures such as coloring and size]p(x ∣ mi)π0(mi)π(mi ∣ x) =∑j p(x ∣ mj)π0(mj)• Problem: (Bayesian) model selection doesn’t

Fin whale<latexit sha1_base64="psq2jvV0jDNdqjwPulyjlu1ZWqw=">AAACAXicbVDLSgNBEOz1GeMr6tHLYBC8GHajoMeAF48RzQOSJcxOZrNjZmaXmVkhLDl59qrf4E28+iV+gn/hJNmDSSxoKKq66e4KEs60cd1vZ2V1bX1js7BV3N7Z3dsvHRw2dZwqQhsk5rFqB1hTziRtGGY4bSeKYhFw2gqGNxO/9USVZrF8MKOE+gIPJAsZwcZK96Ln9Uplt+JOgZaJl5My5Kj3Sj/dfkxSQaUhHGvd8dzE+BlWhhFOx8VuqmmCyRAPaMdSiQXVfjY9dYxOrdJHYaxsSYOm6t+JDAutRyKwnQKbSC96E/E/r5Oa8NrPmExSQyWZLQpTjkyMJn+jPlOUGD6yBBPF7K2IRFhhYmw6c1seo+hcME2q46LNxltMYpk0qxXvolK9uyzXanlKBTiGEzgDD66gBrdQhwYQGMALvMKb8+y8Ox/O56x1xclnjmAOztcvFi2W3Q==</latexit>• Running example: systematics• Goal: learn about evolutionary history of a
• Approach: infer which phylogenetic trees areπ(m1 ∣ x) = 0.8m2

π(m2 ∣ x) = 0.1

Fin whale<latexit sha1_base64="fhdmpccEdfAZWx7deDX0E4BpTnA=">AAACAXicbVDLSgNBEOz1GeMr6tHLYBC8GHajoMeAF48RzQOSJcxOZrNjZmaXmVkhLDl59qrf4E28+iV+gn/hJNmDSSxoKKq66e4KEs60cd1vZ2V1bX1js7BV3N7Z3dsvHRw2dZwqQhsk5rFqB1hTziRtGGY4bSeKYhFw2gqGNxO/9USVZrF8MKOE+gIPJAsZwcZK96JX7ZXKbsWdAi0TLydlyFHvlX66/ZikgkpDONa647mJ8TOsDCOcjovdVNMEkyEe0I6lEguq/Wx66hidWqWPwljZkgZN1b8TGRZaj0RgOwU2kV70JuJ/Xic14bWfMZmkhkoyWxSmHJkYTf5GfaYoMXxkCSaK2VsRibDCxNh05rY8RtG5YJpUx0WbjbeYxDJpViveRaV6d1mu1fKUCnAMJ3AGHlxBDW6hDg0gMIAXeIU359l5dz6cz1nripPPHMEcnK9fF8iW3g==</latexit>

Blue whale
Grey whale
Minke whale.
.
.

always work as we might hope…

3

Reproducibility in model selection[H & Miller 2023]

4

Reproducibility in model selection• Problem: infer phylogeny of
13 whale species frommitochondrial coding DNA[H & Miller 2023]

all
Minke GACCCGAACGTAATAA…ATCCGTTCCCATACTCBlue CACCCCCCCGTACTAT…TGAGTCCGAATTGGAAFin
TGTCTTCTACACTCCA…ACAGGTTGTACGTCAC
Grey GGGTCGCTGTAGACCA…GATACCGCTCTCACAT

4

Reproducibility in model selection• Problem: infer phylogeny of
13 whale species frommitochondrial coding DNA[H & Miller 2023]

all

1st half

Minke GACCCGAACGTAATAA…ATCCGTTCCCATACTCBlue CACCCCCCCGTACTAT…TGAGTCCGAATTGGAAFin
TGTCTTCTACACTCCA…ACAGGTTGTACGTCAC
Grey GGGTCGCTGTAGACCA…GATACCGCTCTCACAT

4

Reproducibility in model selection• Problem: infer phylogeny of
13 whale species frommitochondrial coding DNA[H & Miller 2023]

all

1st half

2nd half

Minke GACCCGAACGTAATAA…ATCCGTTCCCATACTCBlue CACCCCCCCGTACTAT…TGAGTCCGAATTGGAAFin
TGTCTTCTACACTCCA…ACAGGTTGTACGTCAC
Grey GGGTCGCTGTAGACCA…GATACCGCTCTCACAT

4

Reproducibility in model selection• Problem: infer phylogeny of
13 whale species frommitochondrial coding DNA• Compute posterior treeall

1st half

2nd half

Minke GACCCGAACGTAATAA…ATCCGTTCCCATACTCBlue CACCCCCCCGTACTAT…TGAGTCCGAATTGGAAFin
TGTCTTCTACACTCCA…ACAGGTTGTACGTCAC
Grey GGGTCGCTGTAGACCA…GATACCGCTCTCACAT

probabilities based on all,
1st half, and 2nd half

[H & Miller 2023]

4

Reproducibility in model selection• Problem: infer phylogeny of
13 whale species frommitochondrial coding DNA• Compute posterior tree

probabilities based on all,
1st half, and 2nd halfall

1st half

2nd half

Minke GACCCGAACGTAATAA…ATCCGTTCCCATACTCBlue CACCCCCCCGTACTAT…TGAGTCCGAATTGGAAFin
TGTCTTCTACACTCCA…ACAGGTTGTACGTCAC
Grey GGGTCGCTGTAGACCA…GATACCGCTCTCACAT

1st halfall• Compute overlap of 99%
high probability regions[H & Miller 2023]

4

Reproducibility in model selection• Problem: infer phylogeny of
13 whale species frommitochondrial coding DNA• Compute posterior tree

probabilities based on all,
1st half, and 2nd halfall

1st half

2nd half

Minke GACCCGAACGTAATAA…ATCCGTTCCCATACTCBlue CACCCCCCCGTACTAT…TGAGTCCGAATTGGAAFin
TGTCTTCTACACTCCA…ACAGGTTGTACGTCAC
Grey GGGTCGCTGTAGACCA…GATACCGCTCTCACAT

1st halfall• Compute overlap of 99%
high probability regions
all

1st half

overlap

1
0.75
0.5
0.25
0
[H & Miller 2023]

tree 1tree 2tree 3tree 4total

4Reproducibility in model selection• Problem: infer phylogeny of
13 whale species frommitochondrial coding DNA• Compute posterior tree

probabilities based on all,
1st half, and 2nd halfall

1st half

2nd half

Minke GACCCGAACGTAATAA…ATCCGTTCCCATACTCBlue CACCCCCCCGTACTAT…TGAGTCCGAATTGGAAFin
TGTCTTCTACACTCCA…ACAGGTTGTACGTCAC
Grey GGGTCGCTGTAGACCA…GATACCGCTCTCACAT

1st halfall• Compute overlap of 99%
high probability regions
all

1st half

overlap

1
0.75
0.5

0.5

0.25
0
[H & Miller 2023]

tree 1

0.3

0.19

0.0
tree 20.99

tree 3tree 4total

4Reproducibility in model selection• Problem: infer phylogeny of
13 whale species frommitochondrial coding DNA• Compute posterior tree

probabilities based on all,
1st half, and 2nd halfall

1st half

2nd half

Minke GACCCGAACGTAATAA…ATCCGTTCCCATACTCBlue CACCCCCCCGTACTAT…TGAGTCCGAATTGGAAFin
TGTCTTCTACACTCCA…ACAGGTTGTACGTCAC
Grey GGGTCGCTGTAGACCA…GATACCGCTCTCACAT

1st halfall• Compute overlap of 99%
high probability regions
all

1st half

1
0.75
0.5
0
[H & Miller 2023]

0.99 0.99

0.69

0.5

0.25

overlap

0.2
tree 10.0

0.1
tree 2

0.3

0.19

0.0
tree 3tree 4total

4Reproducibility in model selection• Problem: infer phylogeny of
13 whale species frommitochondrial coding DNA• Compute posterior tree

probabilities based on all,
1st half, and 2nd halfall

1st half

2nd half

Minke GACCCGAACGTAATAA…ATCCGTTCCCATACTCBlue CACCCCCCCGTACTAT…TGAGTCCGAATTGGAAFin
TGTCTTCTACACTCCA…ACAGGTTGTACGTCAC
Grey GGGTCGCTGTAGACCA…GATACCGCTCTCACAT

1st halfall• Compute overlap of 99%
high probability regions
all

1st half

1
0.75
0.5
0
[H & Miller 2023]

0.99 0.99

0.69

0.5

0.25

overlap

0.2 0.2
tree 10.0

0.1
tree 2

0.0

0.19

0.19
tree 30.39

0.3
0.0 0.0
tree 4total

4Reproducibility in model selection• Problem: infer phylogeny of
13 whale species frommitochondrial coding DNA• Compute posterior tree

probabilities based on all,
1st half, and 2nd halfall

1st half

2nd half

Minke GACCCGAACGTAATAA…ATCCGTTCCCATACTCBlue CACCCCCCCGTACTAT…TGAGTCCGAATTGGAAFin
TGTCTTCTACACTCCA…ACAGGTTGTACGTCAC
Grey GGGTCGCTGTAGACCA…GATACCGCTCTCACAT

1st halfall• Compute overlap of 99%
high probability regions• 0% overlap = contradiction
all

1st half

1
0.75
0.5
0
[H & Miller 2023]

0.99 0.99

0.69

0.5

0.25

overlap

0.2 0.2
tree 10.0

0.1
tree 2

0.0

0.19

0.19
tree 30.39

0.3
0.0 0.0
tree 4total

4Bayesian phylogenetic inference
may not be reproducible
Problem: infer phylogeny of 13 whale species[H & Miller 2023]

5

Bayesian phylogenetic inference
may not be reproducible
Problem: infer phylogeny of 13 whale species• For some evolutionary models, little to no overlap

Cross-data consistency
all vs 1st half
all vs 2nd half
1st vs 2nd halfevolutionary model[H & Miller 2023]

5

Bayesian phylogenetic inference
may not be reproducible
Problem: infer phylogeny of 13 whale species• For some evolutionary models, little to no overlap

Cross-data consistency
all vs 1st half
all vs 2nd half
1st vs 2nd half

dangerzone
evolutionary model[H & Miller 2023]

5

Bayesian phylogenetic inference
may not be reproducible
Problem: infer phylogeny of 13 whale species• For some evolutionary models, little to no overlap

Cross-data consistency
all vs 1st half
all vs 2nd halfzero overlap

1st vs 2nd half

dangerzone
evolutionary model[H & Miller 2023]

5

Bayesian phylogenetic inference
may not be reproducible
Problem: infer phylogeny of 13 whale species• For some evolutionary models, little to no overlap• (Bayesian) model selection is unstable and not reproducible [Wilcox et al.
2002, Alfaro et al. 2003, Douady et al. 2003, …]

Cross-data consistency
all vs 1st half
all vs 2nd halfzero overlap

1st vs 2nd half

dangerzone
evolutionary model[H & Miller 2023]

5

Bayesian phylogenetic inference
may not be reproducible
Problem: infer phylogeny of 13 whale species• For some evolutionary models, little to no overlap• (Bayesian) model selection is unstable and not reproducible [Wilcox et al.
2002, Alfaro et al. 2003, Douady et al. 2003, …]• Same problem comparing evolutionary models with data fixed
Cross-data consistency

Cross-model consistency
all vs 1st half
all vs 2nd halfzero overlap

1st vs 2nd half

dangerzone
evolutionary model[H & Miller 2023]

evolutionary model
5

Overview

6Overview• Minimal goal: avoid methods that lead to contradictory/non-reproducible inferences

6

Overview• Minimal goal: avoid methods that lead to contradictory/non-reproducible inferences• For example, when…
• getting more data
• slightly tweaking model• using diﬀerent data from same generative process6Overview• Minimal goal: avoid methods that lead to contradictory/non-reproducible inferences• For example, when…
• getting more data
• slightly tweaking model• using diﬀerent data from same generative process• This talk: when and how contradictions can arise in
1. model selection
2. prediction with high-dimensional models
3. unsupervised learning

6

Overview• Minimal goal: avoid methods that lead to contradictory/non-reproducible inferences• For example, when…
• getting more data
• slightly tweaking model• using diﬀerent data from same generative process• This talk: when and how contradictions can arise in
1. model selection
2. prediction with high-dimensional models
3. unsupervised learning• Takeaways:
A. Non-reproducibility can be subtle (and is problem-dependent)
B. Not specific to Bayes
C. Need default, low-cost protective methods that remain statistically eﬃcient
6Part 1: Model Selection[Yang & Zhu 2018, H & Miller 2023]

7

Part 1: Model Selection• Assume data x and models m1 and m2[Yang & Zhu 2018, H & Miller 2023]

7

Part 1: Model Selection• Assume data x and models m1 and m2• Assume they explain the data-generating distribution equally well:𝔼{log p(x ∣ m1)} = 𝔼{log p(x ∣ m2)}

[Yang & Zhu 2018, H & Miller 2023]

7

Part 1: Model Selection• Assume data x and models m1 and m2• Assume they explain the data-generating distribution equally well:𝔼{log p(x ∣ m1)} = 𝔼{log p(x ∣ m2)}

• We’d hope models have equal posterior probability as N → ∞:π(m1 ∣ x) = π(m2 ∣ x) = 1/2[Yang & Zhu 2018, H & Miller 2023]

7

Part 1: Model Selection• Assume data x and models m1 and m2• Assume they explain the data-generating distribution equally well:𝔼{log p(x ∣ m1)} = 𝔼{log p(x ∣ m2)}

• We’d hope models have equal posterior probability as N → ∞:π(m1 ∣ x) = π(m2 ∣ x) = 1/2• However…
Theorem [Yang & Zhu 2019, H & Miller 2023]:d

As N → ∞, π(m1 ∣ x) → Bernoulli(0.5)[Yang & Zhu 2018, H & Miller 2023]

7

Part 1: Model Selection• Assume data x and models m1 and m2• Assume they explain the data-generating distribution equally well:𝔼{log p(x ∣ m1)} = 𝔼{log p(x ∣ m2)}

• We’d hope models have equal posterior probability as N → ∞:π(m1 ∣ x) = π(m2 ∣ x) = 1/2• However…
Theorem [Yang & Zhu 2019, H & Miller 2023]:d

As N → ∞, π(m1 ∣ x) → Bernoulli(0.5)
all posterior mass on asingle, arbitrary model[Yang & Zhu 2018, H & Miller 2023]

7

This is not an abstract concern…

8This is not an abstract concern…

Despite broad agreement that Homo sapiens originated in Africa,
considerable uncertainty surrounds specific models of divergence andmigration across the continent…

8

This is not an abstract concern…

Despite broad agreement that Homo sapiens originated in Africa,
considerable uncertainty surrounds specific models of divergence andmigration across the continent…
…Progress is hampered by a shortage of fossil and genomic data, as wellas variability in previous estimates of divergence times. Here we seek
to discriminate among such models…

8This is not an abstract concern…

Despite broad agreement that Homo sapiens originated in Africa,
considerable uncertainty surrounds specific models of divergence andmigration across the continent…
…Progress is hampered by a shortage of fossil and genomic data, as wellas variability in previous estimates of divergence times. Here we seek
to discriminate among such models…
…We show that model misspecification explains the variation inprevious estimates of divergence time8

A solution: the bagged posteriorx
x*[Douady et al. 2003, Bühlmann 2014, H & Miller 2023]

9

A solution: the bagged posterior• Have data x = (x1, …, xN )x
x*[Douady et al. 2003, Bühlmann 2014, H & Miller 2023]

9

A solution: the bagged posterior• Have data x = (x1, …, xN )• Empirical data distribution ℙN = N−1

N∑n=1 δxnℙN

xx*[Douady et al. 2003, Bühlmann 2014, H & Miller 2023]

9

A solution: the bagged posterior• Have data x = (x1, …, xN )−1

N∑n=1 δxn• Bootstrap dataset x* = (x*1 , …, x*M),xx*• Empirical data distribution ℙN = N
where x*m i.i.d. ∼ ℙN

[Douady et al. 2003, Bühlmann 2014, H & Miller 2023]ℙN

9

A solution: the bagged posterior• Have data x = (x1, …, xN )−1

N∑n=1 δxn• Bootstrap dataset x* = (x*1 , …, x*M),xx*• Empirical data distribution ℙN = N
where x*m i.i.d. ∼ ℙN

[Douady et al. 2003, Bühlmann 2014, H & Miller 2023]ℙN

not always equal to N!

9A solution: the bagged posterior• Have data x = (x1, …, xN )−1

N∑n=1 δxn• Bootstrap dataset x* = (x*1 , …, x*M),xx*• Empirical data distribution ℙN = N
where x*m i.i.d. ∼ ℙN

ℙN

not always equal to N!

• Bagged posterior (a.k.a. BayesBag):
π*(θ ∣ x) = 𝔼{π(θ ∣ x*) ∣ x}

[Douady et al. 2003, Bühlmann 2014, H & Miller 2023]

9

A solution: the bagged posterior• Have data x = (x1, …, xN )−1

N∑n=1 δxn• Bootstrap dataset x* = (x*1 , …, x*M),xx*• Empirical data distribution ℙN = N
where x*m i.i.d. ∼ ℙN

ℙN

not always equal to N!

• Bagged posterior (a.k.a. BayesBag):
π*(θ ∣ x) = 𝔼{π(θ ∣ x*) ∣ x}•Bπ(θ ∣ x*In practice, sample B bootstrap datasets: π*(θ ∣ x) ≈ ∑
(b))
b=1

[Douady et al. 2003, Bühlmann 2014, H & Miller 2023]

9

A solution: the bagged posterior• Have data x = (x1, …, xN )−1

N∑n=1 δxn• Bootstrap dataset x* = (x*1 , …, x*M),xx*• Empirical data distribution ℙN = N
where x*m i.i.d. ∼ ℙN

ℙN

not always equal to N!

• Bagged posterior (a.k.a. BayesBag):
π*(θ ∣ x) = 𝔼{π(θ ∣ x*) ∣ x}•Bπ(θ ∣ x*In practice, sample B bootstrap datasets: π*(θ ∣ x) ≈ ∑
(b))
b=1

‣ Suﬃces to take B = 50 or 100[Douady et al. 2003, Bühlmann 2014, H & Miller 2023]

9

A solution: the bagged posterior• Have data x = (x1, …, xN )−1

N∑n=1 δxn• Bootstrap dataset x* = (x*1 , …, x*M),xx*• Empirical data distribution ℙN = N
where x*m i.i.d. ∼ ℙN

ℙN

not always equal to N!

• Bagged posterior (a.k.a. BayesBag):
π*(θ ∣ x) = 𝔼{π(θ ∣ x*) ∣ x}•Bπ(θ ∣ x*In practice, sample B bootstrap datasets: π*(θ ∣ x) ≈ ∑
(b))
b=1

‣ Suﬃces to take B = 50 or 100‣ Benefits: easy to use, can parallelize across B[Douady et al. 2003, Bühlmann 2014, H & Miller 2023]

9

The bagged posterior is stable• Assume data x and two models m1 and m2• Assume they explain the data-generating distribution equally well:𝔼{log p(x ∣ m1)} = 𝔼{log p(x ∣ m2)}

• We’d hope models have equal poster probability as N → ∞:π(m1 ∣ Y ) = π(m2 ∣ Y ) = 1/2• However…
Theorem [Yang & Zhu 2019, H & Miller 2023]:

all posterior mass on asingle, arbitrary modeld

As N → ∞, π(m1 ∣ x) → Bernoulli(0.5)
Theorem [H & Miller 2023]: As N → ∞d

1. If M = N, then π*(m1 ∣ x) → Uniform(0,1)
P

2. If M/N → 0, then π*(m1 ∣ x) → 1/2[Yang & Zhu 2018, H & Miller 2023]

10

The bagged posterior is stable• Assume data x and two models m1 and m2• Assume they explain the data-generating distribution equally well:𝔼{log p(x ∣ m1)} = 𝔼{log p(x ∣ m2)}

• We’d hope models have equal poster probability as N → ∞:π(m1 ∣ Y ) = π(m2 ∣ Y ) = 1/2• However…
Theorem [Yang & Zhu 2019, H & Miller 2023]:

all posterior mass on asingle, arbitrary modeld

As N → ∞, π(m1 ∣ x) → Bernoulli(0.5)
Theorem [H & Miller 2023]: As N → ∞d

1. If M = N, then π*(m1 ∣ x) → Uniform(0,1)
P

2. If M/N → 0, then π*(m1 ∣ x) → 1/2[Yang & Zhu 2018, H & Miller 2023]

10

The bagged posterior is stable• Assume data x and two models m1 and m2• Assume they explain the data-generating distribution equally well:𝔼{log p(x ∣ m1)} = 𝔼{log p(x ∣ m2)}

• We’d hope models have equal poster probability as N → ∞:π(m1 ∣ Y ) = π(m2 ∣ Y ) = 1/2• However…
Theorem [Yang & Zhu 2019, H & Miller 2023]:

all posterior mass on asingle, arbitrary modeld

As N → ∞, π(m1 ∣ x) → Bernoulli(0.5)
Theorem [H & Miller 2023]: As N → ∞d

1. If M = N, then π*(m1 ∣ x) → Uniform(0,1)
P

2. If M/N → 0, then π*(m1 ∣ x) → 1/2[Yang & Zhu 2018, H & Miller 2023]

10

The bagged posterior is stable• Assume data x and two models m1 and m2• Assume they explain the data-generating distribution equally well:𝔼{log p(x ∣ m1)} = 𝔼{log p(x ∣ m2)}

• We’d hope models have equal poster probability as N → ∞:π(m1 ∣ Y ) = π(m2 ∣ Y ) = 1/2• However…
Theorem [Yang & Zhu 2019, H & Miller 2023]:

all posterior mass on asingle, arbitrary modeld

As N → ∞, π(m1 ∣ x) → Bernoulli(0.5)
Theorem [H & Miller 2023]: As N → ∞d

1. If M = N, then π*(m1 ∣ x) → Uniform(0,1)
P

2. If M/N → 0, then π*(m1 ∣ x) → 1/2[Yang & Zhu 2018, H & Miller 2023]

Recommendation:‣ M = N 0.95 default
0.75
M
=
N
if significant‣
misspecification and/or
many models
10Bagged posterior improves stability
Goal: infer phylogeny of 13 whale species
Cross-data consistency
all vs 1st half
all vs 2nd half
1st vs 2nd half

dangerzone
evolutionary model[H & Miller 2023]

evolutionary model

11

Bagged posterior improves stability
Goal: infer phylogeny of 13 whale species
Cross-data consistency
all vs 1st half
all vs 2nd half
1st vs 2nd half

dangerzone
evolutionary model[H & Miller 2023]

evolutionary model

nonzerooverlap11

Bagged posterior improves stability
Goal: infer phylogeny of 13 whale species
Cross-data consistency
all vs 1st half
all vs 2nd half
1st vs 2nd half

dangerzone
evolutionary model

evolutionary model

nonzerooverlap• Significant overlap between non-MTMAM BayesBag and mixed Bayes• BayesBag dramatically improves cross-model consistency too[H & Miller 2023]

11

Part 2: Prediction with HighDimensional Models

12Part 2: Prediction with HighDimensional Models• Part 1: Bayesian model selection may not be self-consistent (i.e.,
reproducible) if all models are misspecified

12

Part 2: Prediction with HighDimensional Models• Part 1: Bayesian model selection may not be self-consistent (i.e.,
reproducible) if all models are misspecified• But what about parameter estimation?

12

Part 2: Prediction with HighDimensional Models• Part 1: Bayesian model selection may not be self-consistent (i.e.,
reproducible) if all models are misspecified• But what about parameter estimation?
• Usual focus is on pseudo-true parameter (i.e., KL-minimizing parameter)θ ⋆ := arg maxθ 𝔼{log p(x1 ∣ θ)} [Kleijn & van der Vaart 2012, De Blasi & Walker 2012,
Walker 2013, H & Miller 2019, Hoﬀ & Wakefield 2021]

12

Part 2: Prediction with HighDimensional Models• Part 1: Bayesian model selection may not be self-consistent (i.e.,
reproducible) if all models are misspecified• But what about parameter estimation?
• Usual focus is on pseudo-true parameter (i.e., KL-minimizing parameter)θ ⋆ := arg maxθ 𝔼{log p(x1 ∣ θ)} [Kleijn & van der Vaart 2012, De Blasi & Walker 2012,
Walker 2013, H & Miller 2019, Hoﬀ & Wakefield 2021]• This choice is somewhat arbitrary if model misspecified

12

Part 2: Prediction with HighDimensional Models• Part 1: Bayesian model selection may not be self-consistent (i.e.,
reproducible) if all models are misspecified• But what about parameter estimation?
• Usual focus is on pseudo-true parameter (i.e., KL-minimizing parameter)θ ⋆ := arg maxθ 𝔼{log p(x1 ∣ θ)} [Kleijn & van der Vaart 2012, De Blasi & Walker 2012,
Walker 2013, H & Miller 2019, Hoﬀ & Wakefield 2021]• This choice is somewhat arbitrary if model misspecified• We propose an alternative approach (details later if there’s time/interest):

12Part 2: Prediction with HighDimensional Models• Part 1: Bayesian model selection may not be self-consistent (i.e.,
reproducible) if all models are misspecified• But what about parameter estimation?
• Usual focus is on pseudo-true parameter (i.e., KL-minimizing parameter)θ ⋆ := arg maxθ 𝔼{log p(x1 ∣ θ)} [Kleijn & van der Vaart 2012, De Blasi & Walker 2012,
Walker 2013, H & Miller 2019, Hoﬀ & Wakefield 2021]• This choice is somewhat arbitrary if model misspecified• We propose an alternative approach (details later if there’s time/interest):‣ Formalize a generally applicable self-consistency criterion based on
overlap probabilities

12

Part 2: Prediction with HighDimensional Models• Part 1: Bayesian model selection may not be self-consistent (i.e.,
reproducible) if all models are misspecified• But what about parameter estimation?
• Usual focus is on pseudo-true parameter (i.e., KL-minimizing parameter)θ ⋆ := arg maxθ 𝔼{log p(x1 ∣ θ)} [Kleijn & van der Vaart 2012, De Blasi & Walker 2012,
Walker 2013, H & Miller 2019, Hoﬀ & Wakefield 2021]• This choice is somewhat arbitrary if model misspecified• We propose an alternative approach (details later if there’s time/interest):‣ Formalize a generally applicable self-consistency criterion based on
overlap probabilities‣ The posterior can violate criterion (especially in high-dimensional settings)
12

Part 2: Prediction with HighDimensional Models• Part 1: Bayesian model selection may not be self-consistent (i.e.,
reproducible) if all models are misspecified• But what about parameter estimation?
• Usual focus is on pseudo-true parameter (i.e., KL-minimizing parameter)θ ⋆ := arg maxθ 𝔼{log p(x1 ∣ θ)} [Kleijn & van der Vaart 2012, De Blasi & Walker 2012,
Walker 2013, H & Miller 2019, Hoﬀ & Wakefield 2021]• This choice is somewhat arbitrary if model misspecified• We propose an alternative approach (details later if there’s time/interest):‣ Formalize a generally applicable self-consistency criterion based on
overlap probabilities‣ The posterior can violate criterion (especially in high-dimensional settings)‣ The bagged posterior doesn’t violate criterion
12Part 3: Unsupervised Learning ofBona Fide Latent Structure[Box & Hunter 1965, Box 1979, Cox 1990, Lehman 1990, Breiman 2001, Shmueli 2010, Miller & Dunson 2019]13

Part 3: Unsupervised Learning ofBona Fide Latent Structure• Common task: Unsupervised learning of “latent structures”[Box & Hunter 1965, Box 1979, Cox 1990, Lehman 1990, Breiman 2001, Shmueli 2010, Miller & Dunson 2019]13

Part 3: Unsupervised Learning ofBona Fide Latent Structure• Common task: Unsupervised learning of “latent structures”• Phylogenetic inference learns ancestries [Felsenstein 2004][Box & Hunter 1965, Box 1979, Cox 1990, Lehman 1990, Breiman 2001, Shmueli 2010, Miller & Dunson 2019]13

Part 3: Unsupervised Learning ofBona Fide Latent Structure• Common task: Unsupervised learning of “latent structures”• Phylogenetic inference learns ancestries [Felsenstein 2004]• Clustering learns cell types [Lee & McLachlan 2014][Box & Hunter 1965, Box 1979, Cox 1990, Lehman 1990, Breiman 2001, Shmueli 2010, Miller & Dunson 2019]13

Part 3: Unsupervised Learning ofBona Fide Latent Structure• Common task: Unsupervised learning of “latent structures”• Phylogenetic inference learns ancestries [Felsenstein 2004]• Clustering learns cell types [Lee & McLachlan 2014]• Stochastic block model learns communities [Abbe 2018][Box & Hunter 1965, Box 1979, Cox 1990, Lehman 1990, Breiman 2001, Shmueli 2010, Miller & Dunson 2019]13

Part 3: Unsupervised Learning ofBona Fide Latent Structure• Common task: Unsupervised learning of “latent structures”• Phylogenetic inference learns ancestries [Felsenstein 2004]• Clustering learns cell types [Lee & McLachlan 2014]• Stochastic block model learns communities [Abbe 2018]• Latent Dirichlet allocation learns document topics [Blei et al. 2003][Box & Hunter 1965, Box 1979, Cox 1990, Lehman 1990, Breiman 2001, Shmueli 2010, Miller & Dunson 2019]13

Part 3: Unsupervised Learning ofBona Fide Latent Structure• Common task: Unsupervised learning of “latent structures”
explanatory• Phylogenetic inference learns ancestries [Felsenstein 2004]• Clustering learns cell types [Lee & McLachlan 2014]• Stochastic block model learns communities [Abbe 2018]
exploratory• Latent Dirichlet allocation learns document topics [Blei et al. 2003][Box & Hunter 1965, Box 1979, Cox 1990, Lehman 1990, Breiman 2001, Shmueli 2010, Miller & Dunson 2019]13

Part 3: Unsupervised Learning ofBona Fide Latent Structure• Common task: Unsupervised learning of “latent structures”
explanatory• Phylogenetic inference learns ancestries [Felsenstein 2004]• Clustering learns cell types [Lee & McLachlan 2014]• Stochastic block model learns communities [Abbe 2018]
exploratory• Latent Dirichlet allocation learns document topics [Blei et al. 2003][Box & Hunter 1965, Box 1979, Cox 1990, Lehman 1990, Breiman 2001, Shmueli 2010, Miller & Dunson 2019]13

Part 3: Unsupervised Learning ofBona Fide Latent Structure• Common task: Unsupervised learning of “latent structures”
explanatory• Phylogenetic inference learns ancestries [Felsenstein 2004]• Clustering learns cell types [Lee & McLachlan 2014]• Stochastic block model learns communities [Abbe 2018]
exploratory• Latent Dirichlet allocation learns document topics [Blei et al. 2003]• Challenge: how do we recover “true” latent structure when the model is wrong?

[Box & Hunter 1965, Box 1979, Cox 1990, Lehman 1990, Breiman 2001, Shmueli 2010, Miller & Dunson 2019]13

Part 3: Unsupervised Learning ofBona Fide Latent Structure• Common task: Unsupervised learning of “latent structures”
explanatory• Phylogenetic inference learns ancestries [Felsenstein 2004]• Clustering learns cell types [Lee & McLachlan 2014]• Stochastic block model learns communities [Abbe 2018]
exploratory• Latent Dirichlet allocation learns document topics [Blei et al. 2003]• Challenge: how do we recover “true” latent structure when the model is wrong?
• This part: learning bona fide clusters with mixture models[Box & Hunter 1965, Box 1979, Cox 1990, Lehman 1990, Breiman 2001, Shmueli 2010, Miller & Dunson 2019]13

Part 3: Unsupervised Learning ofBona Fide Latent Structure• Common task: Unsupervised learning of “latent structures”
explanatory• Phylogenetic inference learns ancestries [Felsenstein 2004]• Clustering learns cell types [Lee & McLachlan 2014]• Stochastic block model learns communities [Abbe 2018]
exploratory• Latent Dirichlet allocation learns document topics [Blei et al. 2003]• Challenge: how do we recover “true” latent structure when the model is wrong?
• This part: learning bona fide clusters with mixture models• Key ideas: use known causal structure and domain knowledge[Box & Hunter 1965, Box 1979, Cox 1990, Lehman 1990, Breiman 2001, Shmueli 2010, Miller & Dunson 2019]13

Clustering cells by typeflow cytometry data
[cf. Miller & Dunson 2018, Cai, Campbell & Broderick 2021]

14

Clustering cells by type
Standard approach: Gaussian mixture
model with prior on # of components kflow cytometry data
[cf. Miller & Dunson 2018, Cai, Campbell & Broderick 2021]

14

Clustering cells by type
Standard approach: Gaussian mixture
model with prior on # of components kflow cytometry data
[cf. Miller & Dunson 2018, Cai, Campbell & Broderick 2021]

14

Clustering cells by type
Standard approach: Gaussian mixture
model with prior on # of components kflow cytometry data
[cf. Miller & Dunson 2018, Cai, Campbell & Broderick 2021]

More data, more
contradictions!
14Clustering for type discovery

15

Clustering for type discovery• Data-generating process:xn ∼ P∘ =

K∘

∑
k=1

π∘k P∘k‣ Each component is ameaningful type15

Clustering for type discovery• Data-generating process:xn ∼ P∘ =

K∘

∑
k=1

π∘k P∘k‣ Each component is ameaningful type

• Goal: discover these types‣Determine K∘

15

Clustering for type discovery• Data-generating process:xn ∼ P∘ =

K∘

∑
k=1

π∘k P∘kcomponent
distribution‣ Each component is ameaningful type

• Goal: discover these types‣Determine K∘

15

Clustering for type discovery• Data-generating process:xn ∼ P∘ =

K∘

∑
k=1

π∘k P∘kcomponent
distribution‣ Each component is ameaningful type

• Goal: discover these types‣Determine K∘• Assumed mixture model:
xn ∼ Pθ :=

K

∑
k=1

πk Fϕk15

Clustering for type discovery• Data-generating process:xn ∼ P∘ =

K∘

∑
k=1

π∘k P∘kcomponent
distribution‣ Each component is ameaningful type

• Goal: discover these types‣Determine K∘• Assumed mixture model:
xn ∼ Pθ :=

K

∑
k=1

πk Fϕkcomponent
parameter

15

Clustering for type discovery• Data-generating process:xn ∼ P∘ =

K∘

∑
k=1

π∘k P∘kcomponent
distribution‣ Each component is ameaningful type

• Goal: discover these types‣Determine K∘• Assumed mixture model:
xn ∼ Pθ :=

K

∑
k=1

πk Fϕkcomponent
parameter

(K)
θ
=
θ
= (K, π, ϕ1, …, ϕK )Parameters:•

15

Why is type discovery hard?
KπF
k=1 k ϕk

Model: xn ∼ Pθ = ∑K∘
DGP: xn ∼ P∘ = ∑π P
k=1 ∘k ∘k[Miller & Dunson 2018, Cai, Campbell & Broderick 2021]

16

Why is type discovery hard?
KπF
k=1 k ϕk

Model: xn ∼ Pθ = ∑K∘
DGP: xn ∼ P∘ = ∑π P
k=1 ∘k ∘k• Usually, Fϕ ≠ P∘k for all ϕ[Miller & Dunson 2018, Cai, Campbell & Broderick 2021]

16

Why is type discovery hard?
KπF
k=1 k ϕk

Model: xn ∼ Pθ = ∑K∘
DGP: xn ∼ P∘ = ∑π P
k=1 ∘k ∘k• Usually, Fϕ ≠ P∘k for all ϕ• But standard inference
methods ≡ density
estimation[Miller & Dunson 2018, Cai, Campbell & Broderick 2021]

16

Why is type discovery hard?
K
Model: xn ∼ Pθ = ∑πF
k=1 k ϕkN = 100

K∘
DGP: xn ∼ P∘ = ∑π P
k=1 ∘k ∘k• Usually, Fϕ ≠ P∘k for all ϕ• But standard inference
methods ≡ density
estimation• Synthetic example:

skew-normal mixture[Miller & Dunson 2018, Cai, Campbell & Broderick 2021]

16

Why is type discovery hard?
K
Model: xn ∼ Pθ = ∑πF
k=1 k ϕk

N = 500

K∘
DGP: xn ∼ P∘ = ∑π P
k=1 ∘k ∘k• Usually, Fϕ ≠ P∘k for all ϕ• But standard inference
methods ≡ density
estimation• Synthetic example:

skew-normal mixture[Miller & Dunson 2018, Cai, Campbell & Broderick 2021]

16

Why is type discovery hard?
K
Model: xn ∼ Pθ = ∑πF
k=1 k ϕkN = 1000K∘
DGP: xn ∼ P∘ = ∑π P
k=1 ∘k ∘k• Usually, Fϕ ≠ P∘k for all ϕ• But standard inference
methods ≡ density
estimation• Synthetic example:

skew-normal mixture[Miller & Dunson 2018, Cai, Campbell & Broderick 2021]

16

Why is type discovery hard?
K
Model: xn ∼ Pθ = ∑πF
k=1 k ϕk

N = 5000K∘
DGP: xn ∼ P∘ = ∑π P
k=1 ∘k ∘k• Usually, Fϕ ≠ P∘k for all ϕ• But standard inference
methods ≡ density
estimation• Synthetic example:

skew-normal mixture[Miller & Dunson 2018, Cai, Campbell & Broderick 2021]

16

Why is type discovery hard?
K
Model: xn ∼ Pθ = ∑πF
k=1 k ϕkN = 10000

K∘
DGP: xn ∼ P∘ = ∑π P
k=1 ∘k ∘k• Usually, Fϕ ≠ P∘k for all ϕ• But standard inference
methods ≡ density
estimation• Synthetic example:

skew-normal mixture[Miller & Dunson 2018, Cai, Campbell & Broderick 2021]

16

Skew-Gaussian mixture: coarsened posterior[Miller & Dunson 2018, Li & H 2023+]

17

Skew-Gaussian mixture: coarsened posterior
αN+α• One solution: coarsened posterior πα(θ ∣ X) ∝ p(X ∣ θ) π0(θ)[Miller & Dunson 2018, Li & H 2023+]

17

Skew-Gaussian mixture: coarsened posterior
αN+α• One solution: coarsened posterior πα(θ ∣ X) ∝ p(X ∣ θ) π0(θ)
−1
min
KL(P
∣
P
)
≲
α
Roughly,‣θ∘θ[Miller & Dunson 2018, Li & H 2023+]

17

Skew-Gaussian mixture: coarsened posterior
αN+α• One solution: coarsened posterior πα(θ ∣ X) ∝ p(X ∣ θ) π0(θ)
−1
min
KL(P
∣
P
)
≲
α
Roughly,‣θ∘θ[Miller & Dunson 2018, Li & H 2023+]

17

Skew-Gaussian mixture: coarsened posterior
αN+α• One solution: coarsened posterior πα(θ ∣ X) ∝ p(X ∣ θ) π0(θ)
−1
min
KL(P
∣
P
)
≲
α
Roughly,‣θ∘θComputationallycostly!

[Miller & Dunson 2018, Li & H 2023+]

17

Skew-Gaussian mixture: coarsened posterior
αN+α• One solution: coarsened posterior πα(θ ∣ X) ∝ p(X ∣ θ) π0(θ)
−1
min
KL(P
∣
P
)
≲
α
Roughly,‣θ∘θComputationallycostly!

[Miller & Dunson 2018, Li & H 2023+]

17

(Causal) structure matters
Data-generating processK∘
xn ∼ P∘ = ∑k=1 π∘k P∘k18

(Causal) structure matters
Data-generating processK∘
xn ∼ P∘ = ∑k=1 π∘k P∘kStructural causal modelεxx = g(εx)

18(Causal) structure matters
Data-generating process

Structural causal modelK∘
xn ∼ P∘ = ∑k=1 π∘k P∘kεxx = g(εx)“noise” r.v.

18

(Causal) structure matters
Data-generating process

Structural causal modelK∘
xn ∼ P∘ = ∑k=1 π∘k P∘kεxx = g(εx)“noise” r.v.

18

(Causal) structure matters
Data-generating process

Structural causal modelK∘
xn ∼ P∘ = ∑k=1 π∘k P∘kεxx = g(εx)“noise” r.v.

zn ∼ Categorical(π∘)
xn ∣ zn = k ∼ P∘k18

(Causal) structure matters
Data-generating process

Structural causal modelK∘
xn ∼ P∘ = ∑k=1 π∘k P∘kεxx = g(εx)εzz = gz(εz)εxx = gx(z, εx)“noise” r.v.

zn ∼ Categorical(π∘)
xn ∣ zn = k ∼ P∘k18

(Causal) structure matters
Structural causal model

Data-generating processK∘
xn ∼ P∘ = ∑k=1 π∘k P∘kεxx = g(εx)εzz = gz(εz)εxx = gx(z, εx)“noise” r.v.

zn ∼ Categorical(π∘)
xn ∣ zn = k ∼ P∘kAssumed modelzn ∼ Categorical(π)
xn ∣ zn = k ∼ Fϕk
18

(Causal) structure matters
Structural causal model

Data-generating processK∘
xn ∼ P∘ = ∑k=1 π∘k P∘kεxx = g(εx)εzz = gz(εz)εxx = gx(z, εx)“noise” r.v.

zn ∼ Categorical(π∘)
xn ∣ zn = k ∼ P∘kAssumed modelzn ∼ Categorical(π)
xn ∣ zn = k ∼ Fϕk
18

(Causal) structure matters
Structural causal model

Data-generating processK∘
xn ∼ P∘ = ∑k=1 π∘k P∘kεxx = g(εx)εzz = gz(εz)εxx = gx(z, εx)“noise” r.v.

zn ∼ Categorical(π∘)
xn ∣ zn = k ∼ P∘kAssumed modelzn ∼ Categorical(π)
xn ∣ zn = k ∼ Fϕk

Fϕk ≠ P∘k18

Structurally aware model selectionzn ∼ Categorical(π)
xn ∣ zn = k ∼ Fϕk[Li & H 2023+]

KPθ = ∑k=1 πk FϕkK∘
xn ∼ P∘ = ∑k=1 π∘k P∘k19

Structurally aware model selectionzn ∼ Categorical(π)
xn ∣ zn = k ∼ Fϕk

K
Pθ = ∑k=1 πk FϕkK∘
xn ∼ P∘ = ∑k=1 π∘k P∘k−1
Recall:forcoarsened
posterior,assume
min
KL(P
∣
P
)
≲
α
•
θ
∘
θ

[Li & H 2023+]

19

Structurally aware model selectionzn ∼ Categorical(π)
xn ∣ zn = k ∼ Fϕk

K
Pθ = ∑k=1 πk FϕkK∘
xn ∼ P∘ = ∑k=1 π∘k P∘k−1
Recall:forcoarsened
posterior,assume
min
KL(P
∣
P
)
≲
α
•
θ
∘
θ

• But want to coarsen at component level[Li & H 2023+]

19

Structurally aware model selectionzn ∼ Categorical(π)
xn ∣ zn = k ∼ Fϕk

K
Pθ = ∑k=1 πk FϕkK∘
xn ∼ P∘ = ∑k=1 π∘k P∘k−1
Recall:forcoarsened
posterior,assume
min
KL(P
∣
P
)
≲
α
•
θ
∘
θ

• But want to coarsen at component level• Our approach:[Li & H 2023+]

19

Structurally aware model selectionzn ∼ Categorical(π)
xn ∣ zn = k ∼ Fϕk

K
Pθ = ∑k=1 πk FϕkK∘
xn ∼ P∘ = ∑k=1 π∘k P∘k−1
Recall:forcoarsened
posterior,assume
min
KL(P
∣
P
)
≲
α
•
θ
∘
θ

• But want to coarsen at component level• Our approach:‣ Choose divergence 𝒟(P ∣ Q)

[Li & H 2023+]

19

Structurally aware model selectionzn ∼ Categorical(π)
xn ∣ zn = k ∼ Fϕk

K
Pθ = ∑k=1 πk FϕkK∘
xn ∼ P∘ = ∑k=1 π∘k P∘k−1
Recall:forcoarsened
posterior,assume
min
KL(P
∣
P
)
≲
α
•
θ
∘
θ

• But want to coarsen at component level• Our approach:‣ Choose divergence 𝒟(P ∣ Q)
‣ Assumption: minϕ 𝒟(P∘k ∣ Fϕ) ≲ ρ[Li & H 2023+]

19

Structurally aware model selectionzn ∼ Categorical(π)
xn ∣ zn = k ∼ Fϕk

K
Pθ = ∑k=1 πk FϕkK∘
xn ∼ P∘ = ∑k=1 π∘k P∘k−1
Recall:forcoarsened
posterior,assume
min
KL(P
∣
P
)
≲
α
•
θ
∘
θ

• But want to coarsen at component level• Our approach:‣ Choose divergence 𝒟(P ∣ Q)
‣ Assumption: minϕ 𝒟(P∘k ∣ Fϕ) ≲ ρ̂ 1, …, xn ∣ Q) for 𝒟(P ∣ Q)
‣ Need (consistent) estimator 𝒟 (x[Li & H 2023+]

19

Structurally aware model selectionzn ∼ Categorical(π)
xn ∣ zn = k ∼ Fϕk

K
Pθ = ∑k=1 πk FϕkK∘
xn ∼ P∘ = ∑k=1 π∘k P∘k−1
Recall:forcoarsened
posterior,assume
min
KL(P
∣
P
)
≲
α
•
θ
∘
θ

• But want to coarsen at component level• Our approach:‣ Choose divergence 𝒟(P ∣ Q)
‣ Assumption: minϕ 𝒟(P∘k ∣ Fϕ) ≲ ρ̂ 1, …, xn ∣ Q) for 𝒟(P ∣ Q)
‣ Need (consistent) estimator 𝒟 (x‣ Xk(z) := {xn : zn = k} [observations assigned to kth component][Li & H 2023+]

19

Structurally aware model selectionzn ∼ Categorical(π)
xn ∣ zn = k ∼ Fϕk

K
Pθ = ∑k=1 πk FϕkK∘
xn ∼ P∘ = ∑k=1 π∘k P∘k−1
Recall:forcoarsened
posterior,assume
min
KL(P
∣
P
)
≲
α
•
θ
∘
θ

• But want to coarsen at component level• Our approach:‣ Choose divergence 𝒟(P ∣ Q)
‣ Assumption: minϕ 𝒟(P∘k ∣ Fϕ) ≲ ρ̂ 1, …, xn ∣ Q) for 𝒟(P ∣ Q)
‣ Need (consistent) estimator 𝒟 (x‣ Xk(z) := {xn : zn = k} [observations assigned to kth component]‣ Structurally aware loss:K̂ k(z) ∣ Fϕ ) − ρ} + λK
ℒρ(θ ∣ z) = ∑k=1 | Xk(z) | max {0, 𝒟 (Xk

[Li & H 2023+]

19

Structurally aware model selectionzn ∼ Categorical(π)
xn ∣ zn = k ∼ Fϕk

K
Pθ = ∑k=1 πk FϕkK∘
xn ∼ P∘ = ∑k=1 π∘k P∘k−1
Recall:forcoarsened
posterior,assume
min
KL(P
∣
P
)
≲
α
•
θ
∘
θ

• But want to coarsen at component level• Our approach:‣ Choose divergence 𝒟(P ∣ Q)
‣ Assumption: minϕ 𝒟(P∘k ∣ Fϕ) ≲ ρ̂ 1, …, xn ∣ Q) for 𝒟(P ∣ Q)
‣ Need (consistent) estimator 𝒟 (x‣ Xk(z) := {xn : zn = k} [observations assigned to kth component]
AIC-like penalty‣ Structurally aware loss:K̂ k(z) ∣ Fϕ ) − ρ} + λK
ℒρ(θ ∣ z) = ∑k=1 | Xk(z) | max {0, 𝒟 (Xk

[Li & H 2023+]

19

STructurally Aware Robust Estimation(STARE)[Li & H 2023+]

20

STructurally Aware Robust Estimation(STARE)̂ for each K ∈ {1,…, K }1. Estimate model parameters θ(K)max

[Li & H 2023+]

20

STructurally Aware Robust Estimation(STARE)̂ for each K ∈ {1,…, K }1. Estimate model parameters θ(K)max
2. Select ρ[Li & H 2023+]

20

STructurally Aware Robust Estimation(STARE)̂ for each K ∈ {1,…, K }1. Estimate model parameters θ(K)max
2. Select ρ̂ ∣ z (K))
3. Set K̂ρ = arg minK ℒρ(θ(K)[Li & H 2023+]

20

STructurally Aware Robust Estimation(STARE)̂ for each K ∈ {1,…, K }1. Estimate model parameters θ(K)max
2. Select ρ̂ ∣ z (K))
3. Set K̂ρ = arg minK ℒρ(θ(K)• Open question: how to choose ρ[Li & H 2023+]

20

Consistency of STARE[Li & H 2023+]

21

Consistency of STARE(K)
Definition: π⋆(K), ϕ (K), …, ϕ⋆K= asymptotically inferred parameters

given K

[Li & H 2023+]⋆1

21

Consistency of STARE(K)
Definition: π⋆(K), ϕ (K), …, ϕ⋆K= asymptotically inferred parameters

given K⋆1K∘
π P s.t.
Key Assumptions: Exists decomposition of interest P∘ = ∑k=1 ∘k ∘k[Li & H 2023+]

21

Consistency of STARE(K)
Definition: π⋆(K), ϕ (K), …, ϕ⋆K= asymptotically inferred parameters

given K⋆1K∘
π P s.t.
Key Assumptions: Exists decomposition of interest P∘ = ∑k=1 ∘k ∘k(i) model components close when K correct:(K)) < ρKL(P∘k ∣ ∣ Fϕ⋆k

[Li & H 2023+]

21

Consistency of STARE(K)
Definition: π⋆(K), ϕ (K), …, ϕ⋆K= asymptotically inferred parameters

given K⋆1K∘
π P s.t.
Key Assumptions: Exists decomposition of interest P∘ = ∑k=1 ∘k ∘k(i) model components close when K correct:(K)) < ρKL(P∘k ∣ ∣ Fϕ⋆k(ii) smaller mixtures are a poor fit: for π err := ∥π∘ − π⋆(K∘)∥1,
K
(K)
err
(K), P ) > (1 − π
if K < K∘, then dBL( ∑πF)∘
k=1 ⋆k ϕ⋆k[Li & H 2023+]

ρ/2 + π err

21Consistency of STARE(K)
Definition: π⋆(K), ϕ (K), …, ϕ⋆K= asymptotically inferred parameters

given K⋆1K∘
π P s.t.
Key Assumptions: Exists decomposition of interest P∘ = ∑k=1 ∘k ∘k(i) model components close when K correct:(K)) < ρKL(P∘k ∣ ∣ Fϕ⋆k(ii) smaller mixtures are a poor fit: for π err := ∥π∘ − π⋆(K∘)∥1,
K
(K)
err
(K), P ) > (1 − π
if K < K∘, then dBL( ∑πF)∘
k=1 ⋆k ϕ⋆kρ/2 + π err

Theorem [Li & H 2023+]: As N → ∞, Pr(K̂ρ = K∘) → 1.

[Li & H 2023+]

21

Consistency of STARE(K)
Definition: π⋆(K), ϕ (K), …, ϕ⋆K= asymptotically inferred parameters

given K⋆1K∘
π P s.t.
Key Assumptions: Exists decomposition of interest P∘ = ∑k=1 ∘k ∘k(i) model components close when K correct:(K)) < ρKL(P∘k ∣ ∣ Fϕ⋆k(ii) smaller mixtures are a poor fit: for π err := ∥π∘ − π⋆(K∘)∥1,
K
(K)
err
(K), P ) > (1 − π
if K < K∘, then dBL( ∑πF)∘
k=1 ⋆k ϕ⋆kρ/2 + π err

Theorem [Li & H 2023+]: As N → ∞, Pr(K̂ρ = K∘) → 1.

‣ Special case of our general consistency result[Li & H 2023+]

21

Skew-Gaussian mixture: STAREK̂̂ k(z) ∣ F ̂ ) − ρ} + λK
ℒρ(θ ∣ z) = ∑k=1 | Xk(z) | max {0, KL (Xϕk[Li & H 2023+]

22

Skew-Gaussian mixture: STAREK̂̂ k(z) ∣ F ̂ ) − ρ} + λK
ℒρ(θ ∣ z) = ∑k=1 | Xk(z) | max {0, KL (Xϕkpiecewise linear in ρ

[Li & H 2023+]

22

Skew-Gaussian mixture: STAREK̂̂ k(z) ∣ F ̂ ) − ρ} + λK
ℒρ(θ ∣ z) = ∑k=1 | Xk(z) | max {0, KL (Xϕkpiecewise linear in ρ

[Li & H 2023+]

22

Skew-Gaussian mixture: STAREK̂̂ k(z) ∣ F ̂ ) − ρ} + λK
ℒρ(θ ∣ z) = ∑k=1 | Xk(z) | max {0, KL (Xϕkwide, stable region[Li & H 2023+]

22

Skew-Gaussian mixture: STAREK̂̂ k(z) ∣ F ̂ ) − ρ} + λK
ℒρ(θ ∣ z) = ∑k=1 | Xk(z) | max {0, KL (Xϕkwide, stable region[Li & H 2023+]

22

Clustering cells by type: calibrating α

flow cytometry data

[Miller & Dunson 2018]23

Clustering cells by type: calibrating α
α = 200

clustering
accuracy(F measure)[Miller & Dunson 2018]23

Clustering cells by type: calibrating ρ[Li & H 2023+]

24

Clustering cells by type: calibrating ρρ ≈ 1.16
Dataset #clustering
accuracy(F measure)[Li & H 2023+]

24

Clustering cells by type: calibrating ρρ ≈ 1.16
Dataset #clustering
accuracy(F measure)• STARE runtime: 30 min [Python][Li & H 2023+]

24

Clustering cells by type: calibrating ρρ ≈ 1.16
Dataset #clustering
accuracy(F measure)• STARE runtime: 30 min [Python]• Coarsened posterior runtime: 2.5 hours [Julia][Li & H 2023+]

24

A borderline caseDataset 12

ρ = 1.16[Li & H 2023+]

25

Clustering cells by type:
STARE is fast and accurate
Clustering accuracy (F-measure)
Dataset

7

8

9

10

11

12

Runtime

STARE

0.63

0.92

0.94

0.99

0.99

0.98

9 min

Coarsened
posterior

0.67

0.88

0.93

0.99

0.99

0.99

[Miller & Dunson 2018, Li & H 2023+]

26

Closing thoughts27

Closing thoughts• Model misspecification can dramatically aﬀect stability/reproducibility of inferences‣ Ideally, want to default to stable methods that don’t degrade statisticaleﬃciency‣ Examples: Bagged posterior and STARE

27

Closing thoughts• Model misspecification can dramatically aﬀect stability/reproducibility of inferences‣ Ideally, want to default to stable methods that don’t degrade statisticaleﬃciency‣ Examples: Bagged posterior and STARE
ReferencesJ. H. Huggins & J. W. Miller (2023). Reproducible Model Selection Using Bagged
Posteriors. Bayesian Analysis 18(1): 79–104
J. H. Huggins & J. W. Miller (2019). Robust Inference and Model Criticism Using
Bagged Posteriors. arXiv:1912.07104 [stat.ME].
J. H. Huggins & J. W. Miller (2023+). Reproducible Parameter Inference Using Bagged
Posteriors.
J. Li & J. H. Huggins (2023+). Robust, Structurally-Aware Inference with Mixture
Models.
27