Commit 04bdc17e authored by dsbowen's avatar dsbowen
Browse files

Merge branch 'db/update_dirac_delta' into 'master'

Db/update dirac delta

See merge request !2
parents d51c4fba 484688ab
Loading
Loading
Loading
Loading
Loading
+1 −3
Original line number Diff line number Diff line
@@ -38,9 +38,7 @@ typehint:
# Make sphinx docs, run doctests, and serve docs as a web page
.PHONY: docmake
docmake:
	python setup.py build_sphinx\
		--source ${SPHINX_SOURCE_DIR}\
		--build-dir ${SPHINX_BUILD_DIR}
	sphinx-build ${SPHINX_SOURCE_DIR} ${SPHINX_BUILD_DIR}/html
.PHONY: doctest
doctest:
	sphinx-build ${SPHINX_SOURCE_DIR} ${REPORTS_DIR} -b doctest
+0 −1
Original line number Diff line number Diff line
@@ -19,7 +19,6 @@ multiple\_inference.stats
   
      joint_distribution
      mixture
      nonparametric
      quantile_unbiased
      truncnorm
      
+3 −7
Original line number Diff line number Diff line
%% Cell type:markdown id: tags:

# Bayes primer

I designed this notebook to give you a primer on Bayesian analysis: how it works, why you should use it, and how it can change your results. To run Bayesian analysis on your data, click the badge below.

[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gl/dsbowen%2Fconditional-inference/HEAD?urlpath=lab/tree/docs/examples/multiple_inference.ipynb)

First, when should you use Bayesian analysis? You should use Bayesian analysis when comparing 4 or more "things." For example, you should use Bayesian analysis when you run a study comparing the effects of 4 or more treatments or when studying differences between 4 or more groups of people. (The reason we start at 4 instead of 3 or 5 has to do with the [mathematical underpinnings](https://en.wikipedia.org/wiki/James%E2%80%93Stein_estimator) of Bayesian estimators.)

Throughout this notebook, I'll illustrate the importance of Bayesian estimators with an example from [A megastudy of text-based nudges encouraging patients to get vaccinated at an upcoming doctor's appointment](https://www.pnas.org/content/118/20/e2101165118) published in PNAS. The authors partnered with Penn Medicine to send patients one of 19 text messages encouraging them to get a flu vaccine. Using OLS, the authors reported their average text message increased vaccination rates by 2.1 people per hundred compared to the control group. The top-performing message was twice as effective, increasing vaccination rates by a stunning [4.6 people per hundred](https://twitter.com/katy_milkman/status/1362579547401687040).

Many popular media outlets, including the [Economist](https://www.economist.com/by-invitation/2020/11/30/katy-milkman-on-how-to-nudge-people-to-accept-a-covid-19-vaccine), the [Washington Post](https://www.washingtonpost.com/outlook/2021/05/24/nudges-vaccination-psychology-messaging/), [CNBC](https://www.cnbc.com/2021/06/26/return-to-office-and-vaccines-how-companies-can-drum-up-enthusiasm.html), [NPR](https://www.npr.org/2021/05/26/1000616898/the-science-behind-vaccine-incentives), and [CNN](https://kyma.com/cnn-health/2021/06/29/this-simple-text-message-can-encourage-people-to-get-vaccinated-researchers-say/), point to this research as a remarkable example of how behavioral economics can encourage people to get vaccinated and potentially save lives during the COVID-19 pandemic. As [Fortune](https://fortune.com/2021/02/20/covid-vaccine-rollout-getting-people-vaccinated-vaccination-rates-behavioral-nudge-wharto/) reported,

> What they found was eye-opening. Precisely *how* a message was worded had a huge impact on whether the patient ended up getting the shot.

Researchers continue to speculate about why the top-performing message was more successful than the others.

Let's start by looking at the results reported in PNAS.

%% Cell type:code id: tags:

``` python
import warnings
warnings.simplefilter("ignore")

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
import statsmodels.api as sm
from IPython import display
from sklearn.model_selection import RepeatedStratifiedKFold

from multiple_inference.bayes import Improper, Nonparametric, Normal

np.random.seed(123)
sns.set()

display.Image(url="https://www.pnas.org/cms/10.1073/pnas.2101165118/asset/7d1e1f26-cdcd-4d3a-b2a1-167d9d49c6d9/assets/images/large/pnas.2101165118fig01.jpg")
```

%% Cell type:markdown id: tags:

Let's start by downloading the estimated effects.

%% Cell type:code id: tags:

``` python
N_PATIENTS = 47306  # number of participants in the study
CONTROL_VACCINATION_RATE = .42  # vaccination rate in the control condition
XLABEL = "Percentage point increase in flu vaccination"
XLIM = (-.03, .08)

summary_df = pd.read_csv("https://osf.io/download/zqbg2/?view_only=c491df37a33840abbdedda4e60176f34").set_index("condition")
# assume the variance of estimated vaccination rate in the control condition
# is approximately equal to the variance of the estimated vaccination rate
# in the treatment conditions
control_variance = 0.5 * (summary_df.SE ** 2).mean()
# note that the covariance between the estimated treatment effects is the
# variance of the estimated vaccination rate in the control condition
mean, cov = summary_df.Estimate, np.diag(summary_df.SE ** 2 - control_variance) + control_variance
ols_results = Improper(mean, cov, endog_names=XLABEL).fit(title="OLS estimates")
ols_results_plot = ols_results.point_plot()
ols_results_plot.axvline(0, linestyle="--")
ols_results_plot.set_xlim(XLIM)
plt.show()
```

%% Cell type:markdown id: tags:

## How does Bayesian analysis work?

In Bayesian analysis, we start with a prior belief. For example, we might expect that each of the treatments we're about to test will increase vaccination rates by 4 percentage points relative to the control condition. Then we collect data and update our belief. [Bayes' Theorem](https://en.wikipedia.org/wiki/Bayes%27_theorem) is a mathematical formula that tells us how much we should update our prior belief based on the data. The updated belief is called a *posterior*.

Where do we get our prior belief?

Classical Bayes takes the prior as a given. For example, you might have a prior belief based on data from previous studies or a survey of subject matter experts.

However, we can often obtain better estimates by using empirical Bayes to estimate the prior from the data. Estimating the prior using data might sound like a contradiction. By definition, the prior is what you expect *before* seeing any data, so doesn't estimating the prior using data undermine what we're trying to do here?

To understand how [empirical Bayes](https://en.wikipedia.org/wiki/Empirical_Bayes_method) estimates the prior, imagine predicting MLB players' on-base percentage (OBP) next season. We might predict that a player's OBP next season will be the same as his OBP in the previous season. But how can we predict the OBP for a rookie with no batting history? One solution is to predict that the rookie's OBP will be similar to last season's rookies' OBP. In Bayesian terms, we've constructed a prior belief about *next* season's rookies' OBP using data from the *previous* season's rookies' rookies' OBP.

We can apply the same logic to the flu study. Imagine we randomly select one text message and put the data for that treatment in a locked box. What should our prior belief about the effect of this text message be? Empirical Bayes says that our prior belief about the effect of the message we locked in the box should be the average effect of the other 18. We can also use the variability in the effects of the other 18 messages to tell us how confident we should be in our prior, giving us a *prior distribution*.

Empirical Bayes estimators can be parametric or non-parametric. Parametric empirical Bayes assumes the shape of the prior distribution. Nonparametric empirical Bayes does not assume the shape of the prior distribution.

Let's look at the prior from a parametric empirical Bayes estimator assuming a normal prior.

%% Cell type:code id: tags:

``` python
parametric_bayes_model = Normal(mean, cov, endog_names=XLABEL)
prior = parametric_bayes_model.get_marginal_prior(0)
lower, upper = prior.ppf(.025), prior.ppf(.975)
print("Prior 95% CI:", lower, upper)
x = np.linspace(lower, upper)
ax = sns.lineplot(x=x, y=prior.pdf(x))
ax.axvline(prior.mean(), linestyle="--")
ax.set_title("Parametric (normal) empirical Bayes prior")
ax.set_xlabel(XLABEL)
xlim = ax.get_xlim()
plt.show()
```

%% Cell type:markdown id: tags:

According to the parametric empirical Bayes prior, there's a 95% chance that each text message increases vaccination rates by between 0 and 4.1 people per hundred.

Now let's look at the prior from a nonparametric empirical Bayes estimator.

%% Cell type:code id: tags:

``` python
nonparametric_bayes_model = Nonparametric(mean, cov, endog_names=XLABEL)
prior = nonparametric_bayes_model.get_marginal_prior(0)
lower, upper = prior.ppf(.025), prior.ppf(.975)
print("Prior 95% CI:", lower, upper)
ax = sns.lineplot(x=x, y=prior.pdf(x))
ax = sns.lineplot(x=prior.xk, y=prior.pk)
ax.axvline(prior.mean(), linestyle="--")
ax.set_title("Nonparametric empirical Bayes prior")
ax.set_xlabel(XLABEL)
ax.set_xlim(xlim)
plt.show()
```

%% Cell type:markdown id: tags:

According to the nonparametric empirical Bayes prior, there's a 95% chance that each text message increases vaccination rates by between 1.9 and 2.3 people per hundred.

Notice that the nonparametric empirical Bayes prior is narrower than the parametric empirical Bayes prior. This is because the parametric empirical Bayes model accounts for uncertainty in our estimates of the prior parameters. By contrast, nonparametric empirical Bayes priors are often too narrow when estimating only a few treatment effects (in this case, 19). I typically prefer parametric empirical Bayes when estimating fewer than 50 treatment effects.
According to the nonparametric empirical Bayes prior, there's a 95% chance that each text message increases vaccination rates by between 0 and 4.4 people per hundred.

### Summary

In Bayesian analysis, we begin with a prior belief about our treatment effects. We then use data to update our belief according to Bayes' theorem. Our updated belief is called the *posterior*.

The key to good Bayesian analysis is a good prior belief. Classical Bayes takes the prior as a given. Empirical Bayes estimates a prior belief from the data. Parametric empirical Bayes assumes the shape of the prior distribution while nonparametric empirical Bayes does not. Nonparametric empirical Bayes is more flexible than parametric empirical Bayes, but often gives unrealistically narrow confidence intervals when estimating only a few parameters. My rough rule is to use parametric empirical Bayes when estimating fewer than 50 treatment effects.

%% Cell type:markdown id: tags:

## Why should I use Bayesian analysis?

Why use this fancy, complicated Bayesian analysis when you can use standard techniques like OLS? The short answer is that Bayesian estimators make better predictions of the true treatment effects than OLS. We can verify that Bayesian estimators are better using mathematical proofs, out-of-sample testing, and reconstruction plots.

First, let's look at the math. James and Stein (1961) proved that their empirical Bayes estimator *dominates* unbiased estimators like OLS. This means that the James-Stein estimator has a lower expected mean squared error than OLS, regardless of the true treatment effects.

Additionally, unbiased estimators like OLS exaggerate the variability of treatment effects. This fictitious variation makes it seem like treatment effects vary widely even if they do not. [Bayesian estimates "shrink" OLS estimates](https://kiwidamien.github.io/shrinkage-and-empirical-bayes-to-improve-inference.html), meaning that the posterior belief always falls between the OLS estimate and the prior. Bayesian estimators reduce and often eliminate fictitious variation by shrinking the OLS estimates.

A second way to verify that Bayesian estimators make better predictions than OLS is to use [out-of-sample testing](https://en.wikipedia.org/wiki/Cross-validation_(statistics)). To understand out-of-sample testing, imagine we decide to run our experiment twice. After the first experiment, we get both Bayesian and OLS estimates. Then, we use these estimates to predict what will happen in the second experiment. After running the second experiment, we can see which estimator was better.

We may not be able to rerun our experiment, but we can simulate this process by splitting our data in half. We'll use one half of the data (the *training set* or *in-sample data*) to train our models and get Bayesian and OLS estimates. Then, we test how well these estimates matched the other half of the data (the *test set* or *out-of-sample data*).

How can we tell how well our estimates "matched" the test set? We'll measure the mean squared error between the estimated effects on the training set to the OLS treatment effect estimates on the test set. The estimator with the lowest mean squared error is the best.

Below, we repeat this splitting procedure many times to see how our Bayesian estimators stack up against OLS.

%% Cell type:code id: tags:

``` python
# first, we need to approximately reconstruct the dataset
# the dataset is not publicly available
def generate_outcomes(vaccination_rate):
    # create a binary indicator of whether the patient got a vaccine
    n_vaccinations = round(vaccination_rate * n_participants_per_arm)
    return n_vaccinations * [1] + (n_participants_per_arm - n_vaccinations) * [0]

n_participants_per_arm = round(N_PATIENTS / (len(summary_df) + 1))
df = pd.DataFrame()
vaccination_rates = list(summary_df.Estimate + CONTROL_VACCINATION_RATE) + [CONTROL_VACCINATION_RATE]
df["vaccinated"] = np.array([generate_outcomes(vaccination_rate) for vaccination_rate in vaccination_rates]).flatten()
df["treatment"] = np.repeat(list(summary_df.index) + ["control"], n_participants_per_arm)
X = pd.get_dummies(df.drop(columns="vaccinated"), prefix="", prefix_sep="")
X["control"] = 1  # set the control arm as the constant regressor
```

%% Cell type:code id: tags:

``` python
def estimate_mean_and_covariance(index):
    # estimate the OLS means (point estimates) and covariance matrix
    X_subset, y_subset = X.iloc[index], df.vaccinated.iloc[index]
    results = sm.OLS(y_subset, X_subset.astype(int)).fit(cov_type="HC3")
    treatment_coefficients = [i for i, coefficient in enumerate(results.model.exog_names) if coefficient != "control"]
    return (
        results.params[treatment_coefficients],
        results.cov_params().values[treatment_coefficients][:, treatment_coefficients]
    )


def plot_improvement(bayes_mse, title):
    # plot how much the Bayesian model's mean squared error improved upon OLS's test mean squared error
    improvement = ols_mse - bayes_mse
    ax = sns.histplot(x=improvement)
    ax.set_title(title)
    ax.set_xlabel("Reduction in mean squared error compared to OLS")
    return ax


ols_mse, parametric_bayes_mse, nonparametric_bayes_mse = [], [], []
kf = RepeatedStratifiedKFold(n_splits=2, n_repeats=5)
for train_index, test_index in kf.split(df, df.treatment):
    train_mean, train_cov = estimate_mean_and_covariance(train_index)
    test_mean, test_cov = estimate_mean_and_covariance(test_index)
    mean_squared_error = lambda model_cls: (
        (model_cls(train_mean, train_cov).fit().params - test_mean) ** 2
    ).mean()

    ols_mse.append(mean_squared_error(Improper))
    parametric_bayes_mse.append(mean_squared_error(Normal))
    nonparametric_bayes_mse.append(mean_squared_error(Nonparametric))

ols_mse, parametric_bayes_mse, nonparametric_bayes_mse = (
    np.array(ols_mse), np.array(parametric_bayes_mse), np.array(nonparametric_bayes_mse)
)
plot_improvement(parametric_bayes_mse, "Parametric empirical Bayes vs. OLS")
plt.show()

plot_improvement(nonparametric_bayes_mse, "Nonparametric empirical Bayes vs. OLS")
plt.show()

pd.DataFrame({
    "Model": ["OLS", "Parametric empirical Bayes", "Nonparametric empirical Bayes"],
    "Mean squared error": [
        mse.mean() for mse in (ols_mse, parametric_bayes_mse, nonparametric_bayes_mse)
    ]
})
```

%% Cell type:markdown id: tags:

Our out-of-sample analysis suggests that Bayesian estimators outperform OLS.

The third way to verify that Bayesian estimators are better than OLS is to look at reconstruction plots. Reconstruction plots are the most intuitive demonstration that Bayesian estimators are superior, although out-of-sample testing is more rigorous.

Reconstruction plots answer the following question: If these estimates are correct and we reran our experiment, how similar would the distribution of estimates in the second experiment be to the distribution of estimates in the original experiment (i.e., the experiment we actually ran)? Ideally, the distribution of estimates we would expect to see if we reran the experiment should match the distribution of estimates we saw in the original.

How do we get the distribution of estimates we would expect to see if we reran the experiment? Unfortunately, this question takes us into Ph.D.-level statistics territory, so I'll refer curious and ambitious readers to the [Wikipedia entry on Gibbs Sampling](https://en.wikipedia.org/wiki/Gibbs_sampling) for more detail.

Below are reconstruction plots for the OLS and Bayesian estimators. The original estimates are the orange x's. The distribution of estimates we would expect to see if we reran the experiment is in blue. Ideally, the blue dots should overlap with the orange x's.

%% Cell type:code id: tags:

``` python
def make_reconstruction_plot(results, title=None, xlim=None):
    ax = results.reconstruction_point_plot(title=title)
    if xlim:
        ax.set_xlim(xlim)
    ax.axvline(0, linestyle="--")
    return ax


ols_reconstruction_plot = make_reconstruction_plot(ols_results, "OLS reconstruction plot")
xlim = ols_reconstruction_plot.get_xlim()
plt.show()

parametric_results = parametric_bayes_model.fit()
make_reconstruction_plot(parametric_results, title="Parametric empirical Bayes reconstruction plot", xlim=xlim)
plt.show()

nonparametric_results = nonparametric_bayes_model.fit()
make_reconstruction_plot(nonparametric_results, title="Nonparametric empirical Bayes reconstruction plot", xlim=xlim)
plt.show()
```

%% Cell type:markdown id: tags:

Notice the blue dots are more spread out than the orange x's in the OLS reconstruction plot. This confirms that OLS suffers from fictitious variation. By contrast, the blue dots are on top of the orange x's in the Bayesian reconstruction plots. This confirms that Bayesian estimators reliably estimate the distribution of treatment effects.

### Summary

Bayesian analysis is better than traditional techniques like OLS because it makes more accurate predictions of the true treatment effects. We verified this using mathematical proofs, out-of-sample testing, and reconstruction plots.

%% Cell type:markdown id: tags:

## How much can Bayesian analysis change my results?

Maybe you're thinking, "Okay, I'm convinced that Bayesian models are better than OLS, but how different are they? Maybe they'll shrink the OLS estimates slightly but is the difference significant? Can Bayesian estimators fundamentally change our understanding of scientific research?"

To understand the impact of Bayesian analysis, let's see how OLS and Bayesian estimates compare for our flu study. As a reminder, the common perception of this study's results, both in popular media and in academic circles, is that the ability of a text to increase vaccination rates critically depends on its phrasing. Fortune best sums up this perception.

> What they found was eye-opening. Precisely *how* a message was worded had a huge impact on whether the patient ended up getting the shot.

Now that we've verified that Bayesian estimators outperform OLS, let's plot the OLS and Bayesian results.

%% Cell type:code id: tags:

``` python
def point_plot(results, **kwargs):
    ax = results.point_plot(**kwargs)
    ax.set_xlim(XLIM)
    ax.axvline(0, linestyle="--")
    return ax


point_plot(ols_results, title="OLS estimates")
plt.show()

point_plot(parametric_results, title="Parametric empirical Bayes estimates")
plt.show()

point_plot(nonparametric_results, title="Nonparametric empirical Bayes estimates")
plt.show()
```

%% Cell type:markdown id: tags:

Looking at the OLS plot, we get the impression that we've identified certain text messages that outperform others. The Bayesian plots show us that this perception is incorrect. According to the Bayesian models, the treatment effects are extremely similar.

Side note: Remember how we saw earlier that the nonparametric empirical Bayes prior was unrealistically narrow? The narrow prior leads to a narrow posterior (see the nonparametric empirical Bayes plot above). This is why parametric empirical Bayes is often better when analyzing only a few treatments.

According to OLS, the top-performing text message increases vaccination rates by 25 people per thousand compared to the average message. According to parametric empirical Bayes, the top-performing message increases vaccination rates by only four people per thousand. Nonparametric empirical Bayes suggests that the top-performing message increases vaccination rates by less than one person per thousand.
According to OLS, the top-performing text message increases vaccination rates by 25 people per thousand compared to the average message. According to parametric empirical Bayes, the top-performing message increases vaccination rates by only four people per thousand. Nonparametric empirical Bayes suggests that the top-performing message increases vaccination rates by nine people per thousand.

%% Cell type:code id: tags:

``` python
print(
    "OLS: Increase in vaccination rates using the top performing treatment versus the average treatment:",
    1000 * (ols_results.params[0] - ols_results.params.mean()),
    "per 1,000"
)
print(
    "Parametric empirical Bayes: Increase in vaccination rates using the top performing treatment versus the average treatment:",
    1000 * (parametric_results.params[0] - parametric_results.params.mean()),
    "per 1,000"
)
print(
    "Nonparametric empirical Bayes: Increase in vaccination rates using the top performing treatment versus the average treatment:",
    1000 * (nonparametric_results.params[0] - nonparametric_results.params.mean()),
    "per 1,000"
)
```

%% Cell type:markdown id: tags:

According to OLS, the phrasing matters tremendously. Our Bayesian models again show us that this picture is incorrect.

In sum, texting patients a reminder to get a flu vaccine boosts vaccination rates by about 2.1 people per hundred. Beyond the mere act of texting a reminder, there's no evidence that the phrasing of the text messages used in this study has a practically significant effect on vaccination rates.

## Conclusion

Bayesian analysis can significantly change our understanding of scientific research. We illustrated this by re-analyzing data from a highly-regarded study. Using Bayesian analysis, we showed that the study's original conclusion vastly overstated the effect of the top-performing treatment compared to the average treatment.

Congratulations for sticking with this primer until the end! We've explained how Bayesian analysis works, why you should use it, and how it can impact your results.

%% Cell type:code id: tags:

``` python
```
+169 −123
Original line number Diff line number Diff line
@@ -18,36 +18,105 @@ Notes:
"""
from __future__ import annotations

from itertools import product
from typing import Any
from typing import Any, Union

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import minimize_scalar
from scipy.stats import loguniform, norm, rv_continuous
from sklearn.cluster import KMeans
from sklearn.model_selection import check_cv
from sklearn.neighbors import KernelDensity
import seaborn as sns
from scipy.stats import norm, rv_discrete

from .base import BayesBase, BayesResults, ColumnType

# x-values of the Dirac-Delta prior extend some number of standard deviations from the
# min and max values of the conventionally-estimated means
STDS_FROM_MIN_MAX = 2
# minimum number of parameters of the Dirac-Delta prior when ``num_fit`` is "auto"
MIN_PARAMS = 5
# base parameter used for inferring the number of parameters when ``num_fit`` is "auto"
LOG_BASE = 1.3
# epsilon to add to the posterior rvs to break ties (because the posteriors are discrete)
EPS = 1e-8


class NonparametricResults(BayesResults):
    @property
    def _posterior_rvs(self):
        if hasattr(self, "_cached_posterior_rvs"):
            return self._cached_posterior_rvs

        self._cached_posterior_rvs = super()._posterior_rvs
        self._cached_posterior_rvs + EPS * np.random.rand(
            *self._cached_posterior_rvs.shape
        )
        return self._cached_posterior_rvs

    def line_plot(
        self,
        column: ColumnType = None,
        alpha: float = 0.05,
        title: str = None,
        yname: str = None,
        ax=None,
    ):
        """Create a line plot of the prior, conventional, and posterior estimates.

        Args:
            column (ColumnType, optional): Selected parameter. Defaults to None.
            alpha (float, optional): Sets the plot width. 0 is as wide as possible, 1 is
                as narrow as possible. Defaults to .05.
            title (str, optional): Plot title. Defaults to None.
            yname (str, optional): Name of the dependent variable. Defaults to None.
            ax (AxesSubplot, optional): Axis to write on.

        Returns:
            AxesSubplot: Plot.
        """
        index = self.model.get_index(column)
        prior = self.model.get_marginal_prior(index)
        posterior = self.marginal_distributions[index]
        conventional = norm(
            self.model.mean[index], np.sqrt(self.model.cov[index, index])
        )
        xlim = np.array(
            [
                dist.ppf([alpha / 2, 1 - alpha / 2])
                for dist in (prior, conventional, posterior)
            ]
        ).T
        mask = (xlim[0].min() <= prior.xk) & (prior.xk <= xlim[1].max())
        palette = sns.color_palette()
        if ax is None:
            _, ax = plt.subplots()

        sns.lineplot(x=prior.xk[mask], y=prior.pk[mask], label="prior", ax=ax)
        ax.axvline(prior.mean(), linestyle="--", color=palette[0])

        conventional_pmf = conventional.pdf(prior.xk[mask])
        conventional_pmf /= conventional_pmf.sum()
        sns.lineplot(x=prior.xk[mask], y=conventional_pmf, label="conventional")
        ax.axvline(conventional.mean(), linestyle="--", color=palette[1])

        sns.lineplot(
            x=posterior.xk[mask], y=posterior.pk[mask], label="posterior", ax=ax
        )
        ax.axvline(posterior.mean(), linestyle="--", color=palette[2])

        ax.set_title(title or self.model.exog_names[index])
        ax.set_xlabel(yname or self.model.endog_names)

from ..stats import mixture, nonparametric
from .base import BayesBase
        return ax


class Nonparametric(BayesBase):
    """Bayesian model with a nonparametric Dirac delta prior.

    Args:
        num (int, optional): Number of parameters to fit for the prior. Defaults to 100.
        n_clusters (int, optional): Number of clusters to use for featurized
            estimation. Defaults to 1.
        cv (int, optional): Determines the cross validation splitting strategy (input to
            ``sklearn.model_selection.check_cv``). Defaults to 5.
        rtol (float, optional): Relative tolerance stopping criteria for expectation
            maximization. The EM algorithm terminates when the relative improvement
            between iterations falls below this threshold. Defaults to .99.
        max_iter (int, optional): Maximum number of EM iterations. Defaults to 100.
        bandwidth_rvs_size (int, optional): Number of bandwidth values to try when
            tuning the kernel density estimator in between EM iterations to smooth the
            prior. Defaults to 32.
        num_fit (Union[str, int], optional): Number of PMF parameters to fit in the Dirac-Delta prior. Defaults to "auto".
        num_interp (int, optional): Number of PMF parameters to interpolate for the Dirac-Delta prior. Defaults to 1_000.
        lr (float, optional): Learning rate for the Adam optimizer. Defaults to 1e-3.
        betas (tuple[float, float], optional): Beta parameters for Adam. Defaults to (0.9, 0.99).
        eps (float, optional): Epsilon parameter for Adam. Defaults to 1e-8.
        train_iter (int, optional): Number of iterations of gradient descent. Defaults to 1_000.

    Examples:

@@ -69,124 +138,101 @@ class Nonparametric(BayesBase):
            =======================================
                coef pvalue (1-sided) [0.025 0.975]
            ---------------------------------------
            x0 0.837            0.058 -0.197  1.962
            x1 1.200            0.017  0.080  2.968
            x2 1.873            0.003  0.410  3.978
            x3 3.021            0.000  0.927  4.969
            x4 4.037            0.000  1.947  5.858
            x5 4.938            0.000  3.148  7.077
            x6 6.026            0.000  3.981  8.056
            x7 7.092            0.000  5.030  8.573
            x8 7.765            0.000  6.097  8.912
            x9 8.171            0.000  6.942  9.193
            x0 0.535            0.284 -1.068  2.295
            x1 1.340            0.073 -0.433  3.136
            x2 2.191            0.009  0.367  4.039
            x3 3.073            0.001  1.229  4.982
            x4 4.013            0.000  2.111  5.946
            x5 4.987            0.000  3.054  6.889
            x6 5.927            0.000  4.018  7.771
            x7 6.809            0.000  4.961  8.633
            x8 7.660            0.000  5.864  9.433
            x9 8.465            0.000  6.705 10.068
            ===============
            Dep. Variable y
            ---------------
    """

    _results_cls = NonparametricResults

    def __init__(
        self,
        *args: Any,
        num: int = 100,
        n_clusters: int = 1,
        cv=5,
        rtol: float = 0.99,
        max_iter: int = 100,
        bandwidth_rvs_size: int = 32,
        num_fit: Union[str, int] = "auto",
        num_interp: int = 1_000,
        lr: float = 1e-3,
        betas: tuple[float, float] = (0.9, 0.99),
        eps: float = 1e-8,
        train_iter: int = 1_000,
        **kwargs: Any,
    ):
        super().__init__(*args, **kwargs)
        std = self.mean.std()
        lower, upper = self.mean.min() - 2 * std, self.mean.max() + 2 * std
        # (num,) array of values over which the prior is defined
        self._values = np.linspace(lower, upper, num)
        # (num, n_clusters) probability mass function
        self._pmf_values = np.full((num, n_clusters), 1 / num)
        # (# params, n_clusters) mixture weights for each parameter
        self._mixture_weights = KMeans(n_clusters).fit_transform(self.X)
        if (self._mixture_weights == 0).all():
            self._mixture_weights = np.ones(self._mixture_weights.shape)
        self._mixture_weights = (
            self._mixture_weights.T / self._mixture_weights.sum(axis=1)
        ).T

        def loss(value, index, cluster):
            factor = (1 - value) / (1 - self._pmf_values[index, cluster])
            self._pmf_values[:, cluster] *= factor
            self._pmf_values[index, cluster] = value
            arr = self._mixture_weights * (conditional_pdf @ self._pmf_values)
            return -np.log(arr.sum(axis=1)).sum()

        # density function of the conventional estimates evaluated at self._values
        conditional_pdf = [
            norm.pdf(self._values, mean_i, np.sqrt(variance_i))
            for mean_i, variance_i in zip(self.mean, self.cov.diagonal())
        ]
        conditional_pdf = np.array(conditional_pdf)
        # fit the prior using an EM algorithm
        prev_loss, current_loss, i = np.inf, None, 0
        values = self._values.reshape(-1, 1)
        index_cluster = list(product(np.arange(num), np.arange(n_clusters)))
        index_cluster = np.array(index_cluster).astype(int)
        cv = check_cv(cv)
        cv.shuffle = True
        for i in range(max_iter):
            # optimize each value of ``self._pmf_values``
            np.random.shuffle(index_cluster)
            for index, cluster in index_cluster:
                current_loss = minimize_scalar(
                    loss, bounds=(0, 1), method="bounded", args=(index, cluster)
                ).fun

            # smooth the PMF using a kernel density estimator
            cv.random_state = i
            for cluster in range(n_clusters):
                pmf_values = self._pmf_values[:, cluster]
                mean = np.average(self._values, weights=pmf_values)
                std = np.sqrt(
                    np.average((self._values - mean) ** 2, weights=pmf_values)
                )
                bandwidth_rvs = loguniform(0.1 * std, 2 * std).rvs(bandwidth_rvs_size)
                best_score = -np.inf
                for bandwidth in bandwidth_rvs:
                    for train_index, test_index in cv.split(values):
                        X_train, X_test = values[train_index], values[test_index]
                        weight_train = pmf_values[train_index]
                        weight_test = pmf_values[test_index]
                        weight_train /= weight_train.sum()
                        weight_test /= weight_test.sum()
                        kde = KernelDensity(bandwidth=bandwidth).fit(
                            X_train, sample_weight=weight_train
        if num_fit == "auto":
            # this is a very simple heuristic that seems to work well for many datasets
            # ideally, ``num_fit`` should be tuned
            std_ratio = std / np.sqrt(self.cov.diagonal()).mean()
            num_fit = max(
                MIN_PARAMS,
                int(
                    np.log(self.n_params) / LOG_BASE / (max(0, np.log10(std_ratio)) + 1)
                ),
            )
                        score = (weight_test * kde.score_samples(X_test)).mean()
                        if score > best_score:
                            best_score, best_bandwidth = score, bandwidth

                kde = KernelDensity(bandwidth=best_bandwidth).fit(
                    values, sample_weight=pmf_values
        # x-values of the Dirac-Delta prior
        x = np.linspace(
            self.mean.min() - STDS_FROM_MIN_MAX * std,
            self.mean.max() + STDS_FROM_MIN_MAX * std,
            num=num_fit,
        )
                self._pmf_values[:, cluster] = np.exp(kde.score_samples(values))

            self._pmf_values /= self._pmf_values.sum(axis=0)
            if current_loss / prev_loss > rtol:
                break
            prev_loss = current_loss

        # fit a nonparametric distribution for each cluster
        self._cluster_distributions = [
            nonparametric((self._values, self._pmf_values[:, i]))
            for i in range(n_clusters)
        # compute (n_means, num) array f(m | \mu)
        scale = np.sqrt(self.cov.diagonal())
        # note that norm.logpdf(m, loc=mu) = norm.logpdf(mu, loc=m) because the normal
        # distribution is symmetric
        log_pdf_mean_given_mu = np.array(
            [
                norm.logpdf(x, loc=loc_i, scale=scale_i)
                for loc_i, scale_i in zip(self.mean, scale)
            ]
        )
        pdf_mean_given_mu = np.exp(log_pdf_mean_given_mu - log_pdf_mean_given_mu.max())

        # optimize the log(pmf) values of the Dirac-Delta prior to maximize f(m) using Adam
        # optimizer
        pmf = pdf_mean_given_mu.mean(axis=0)
        # initialize log_pmf to the optimal values if the means were estimated without noise
        log_pmf = np.log(pmf / pmf.sum())
        momentum, velocity = 0, 0
        for _ in range(train_iter):
            pmf = np.exp(log_pmf)
            prod = (pmf * pdf_mean_given_mu).sum(axis=1, keepdims=True)
            grad = pmf / (1 - pmf) * (1 - (pdf_mean_given_mu / prod).mean(axis=0))
            momentum = betas[0] * momentum + (1 - betas[0]) * grad
            velocity = betas[1] * velocity + (1 - betas[1]) * grad**2
            momentum_hat = momentum / (1 - betas[0])
            velocity_hat = velocity / (1 - betas[1])
            log_pmf -= lr * momentum_hat / (np.sqrt(velocity_hat) + eps)
            pmf = np.exp(log_pmf)
            log_pmf = np.log(pmf / pmf.sum())

        # interpolate between the fitted values to get a smoother prior
        self._values = np.linspace(x[0], x[-1], num=num_interp)
        pmf = np.interp(self._values, x, pmf)
        self._prior = rv_discrete(values=(self._values, pmf / pmf.sum()))

        # cache values for the computing the posterior
        self._logpdf_mean_given_mu = np.array(
            [
                norm.logpdf(self._prior.xk, loc=loc_i, scale=scale_i)
                for loc_i, scale_i in zip(self.mean, scale)
            ]
        )

    def _get_marginal_prior(self, index: int) -> rv_continuous:
        if len(self._cluster_distributions) == 1:
            return self._cluster_distributions[0]

        return mixture(self._cluster_distributions, self._mixture_weights[index])
    def _get_marginal_prior(self, index: int) -> rv_discrete:
        return self._prior

    def _get_marginal_distribution(self, index: int) -> rv_continuous:
        pmf = (self._pmf_values * self._mixture_weights[index]).sum(axis=1)
        logpmf = np.log(pmf) + norm.logpdf(
            self._values, self.mean[index], np.sqrt(self.cov[index, index])
        )
        return nonparametric((self._values, np.exp(logpmf - logpmf.max())))
    def _get_marginal_distribution(self, index: int) -> rv_discrete:
        log_pmf = self._logpdf_mean_given_mu[index] + np.log(self._prior.pk)
        pmf = np.exp(log_pmf - log_pmf.max())
        return rv_discrete(values=(self._prior.xk, pmf / pmf.sum()))
+2 −144
Original line number Diff line number Diff line
@@ -3,18 +3,14 @@
from __future__ import annotations

import warnings
from typing import Any, Callable, List, Sequence, Tuple, Union
from typing import Any, List, Sequence, Tuple, Union

import numpy as np
from scipy.integrate import quad
from scipy.interpolate import CubicSpline, PchipInterpolator
from scipy.misc import derivative
from scipy.optimize import NonlinearConstraint, fsolve, minimize, minimize_scalar
from scipy.optimize import NonlinearConstraint, fsolve, minimize
from scipy.stats import norm, rv_continuous, truncnorm as truncnorm_base


from .base import Numeric1DArray
from .utils import weighted_quantile


class joint_distribution:
@@ -118,144 +114,6 @@ class mixture(rv_continuous):
        return np.sqrt(self.var())


class nonparametric(rv_continuous):
    """Nonparametric distribution.

    Args:
        values (tuple[np.array, np.array]): (n,) array of x values, (n,) array of the
            probability mass function evaluated at x.
        moment_approximation_samples (int): Number of samples to use when numerically
            approximating moments.

    Attributes:
        xk (np.ndarray): (n,) array of x values.
        pk (np.ndarray): (n,) array of the probability mass function evaluated at x.

    Notes:
        This distribution interpolates between the probability mass function to
        "continuize" the discrete function.
    """

    def __init__(self, values, moment_approximation_samples: int = 50, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.xk, self.pk = np.array(values[0], float), np.array(values[1], float)
        self.moment_approximation_samples = moment_approximation_samples
        if len(self.xk) > 3:
            self._cubic_spline = CubicSpline(self.xk, self.pk)
            self._pchip_interpolater = PchipInterpolator(self.xk, self.pk)

            # put as little weight as possible on the PCHIP interpolater such that the pdf
            # is always above 0
            self._pchip_weight = 0
            min_pdf_result = minimize_scalar(
                self._interpolate, bounds=(self.xk[0], self.xk[-1])
            )
            while min_pdf_result.fun < 0:
                # pdf goes below 0, need to increase weight on the PCHIP interpolater
                pchip_min_fun = self._pchip_interpolater(min_pdf_result.x)
                pchip_weight = min_pdf_result.x / (min_pdf_result.x - pchip_min_fun)
                self._pchip_weight += pchip_weight * (1 - self._pchip_weight)
                min_pdf_result = minimize_scalar(
                    self._interpolate, bounds=(self.xk[0], self.xk[-1])
                )
        else:
            self._cubic_spline = None
            self._pchip_interpolater = None
            self._pchip_weight = None

        # determine the scale so the PDF integrates to 1
        self._scale = 1
        self._scale = 1 / quad(self._pdf, self.xk[0], self.xk[-1])[0]

        # cache evaluations of the CDF
        # this vastly speeds up the ppf and rvs methods
        self._cdf_x = np.linspace(self.xk[0], self.xk[-1], num=10000)
        pdf = self._pdf(self._cdf_x)
        self._cdf_cache = pdf.cumsum()
        self._cdf_cache /= self._cdf_cache[-1]

        super().__init__(*args, **kwargs)

    def _interpolate(self, x: np.ndarray) -> np.ndarray:
        """Apply interpolation to get the PDF evaluated at x.

        Args:
            x (np.ndarray): Points at which to evaluate the PDF.

        Returns:
            np.ndarray: PDF evaluated at x.
        """
        if len(self.xk) < 4:
            return np.interp(x, self.xk, self.pk)

        if self._pchip_weight == 0:
            return self._cubic_spline(x)

        return self._pchip_weight * self._pchip_interpolater(x) + (
            1 - self._pchip_weight
        ) * self._cubic_spline(x)

    def _pdf(self, x: np.ndarray) -> np.ndarray:
        x = np.atleast_1d(x)
        pdf = np.zeros(len(x))
        in_range = (self.xk[0] < x) & (x < self.xk[-1])
        pdf[in_range] = self._interpolate(x[in_range])
        return self._scale * pdf

    def _cdf(self, x: np.ndarray) -> np.ndarray:
        x = np.atleast_1d(x)
        cdf = np.zeros(len(x))
        cdf[x >= self.xk[-1]] = 1
        in_range = (self.xk[0] < x) & (x < self.xk[-1])
        cdf[in_range] = np.interp(x[in_range], self._cdf_x, self._cdf_cache)
        return cdf

    def ppf(self, x: np.ndarray) -> np.ndarray:
        x = np.atleast_1d(x)
        ppf = super()._ppf(x)
        ppf[ppf == np.inf] = self.xk[-1]
        ppf[ppf == -np.inf] = self.xk[0]
        return ppf[0] if len(x) == 1 else ppf

    def moment(self, func: Callable[[np.ndarray], np.ndarray]) -> float:
        """Compute a moment.

        Args:
            func (Callable[[np.ndarray], np.ndarray]): Moment function that takes
                ``self.xk`` and returns an array of the same shape.

        Returns:
            float: Moment.
        """
        x = np.linspace(self.xk[0], self.xk[-1], num=self.moment_approximation_samples)
        return sum(self.pdf(x) * func(x) * (x[1] - x[0]))

    def mean(self) -> float:
        """Compute the mean.

        Returns:
            float: Mean.
        """
        return self.moment(lambda x: x)

    def var(self) -> float:
        """Compute the variance.

        Returns:
            float: Variance.
        """
        mean = self.mean()
        return self.moment(lambda x: (x - mean) ** 2)

    def std(self) -> float:
        """Compute the standard deviation.

        Returns:
            float: Standard deviation.
        """
        return np.sqrt(self.var())


class quantile_unbiased(rv_continuous):  # pylint: disable=invalid-name
    """Conditional quantile-unbiased distribution.

Loading