\subsubsection*{\bibname}}
\bibliographystyle{unsrtnat}

\usepackage{mathtools}
\usepackage{booktabs}
\usepackage{tikz}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{hyperref}
\usepackage{url}
\usepackage{booktabs}
\usepackage{amsfonts}
\usepackage{nicefrac}
\usepackage{microtype}
\usepackage{amsmath}
\usepackage{subcaption}
\usepackage{enumitem}
\usepackage{commath}

\usepackage{caption}
\usepackage{algorithm}
\usepackage[noend]{algpseudocode}
\usepackage{wrapfig}
\usepackage{enumitem}
\usepackage{xr}

\usepackage{amsthm}
\usepackage{thmtools}
\usepackage{thm-restate}

\newtheorem{theorem}{Theorem}
\newtheorem{lemma}{Lemma}
\newtheorem{definition}{Definition}
\newtheorem{remark}{Remark}
\newtheorem{condition}{Condition}
\newtheorem{claim}{Claim}
\usepackage[normalem]{ulem}

\usepackage{tikz}
\usetikzlibrary{bayesnet}

\DeclareMathOperator*{\argmin}{argmin}
\DeclareMathOperator*{\argmax}{argmax}
\DeclareMathOperator*{\Var}{Var}
\DeclareMathOperator*{\Cov}{Cov}
\DeclareMathOperator*{\Lap}{Lap}
\DeclareMathOperator*{\priv}{priv}
\DeclareMathOperator*{\thetapriv}{\hat{\theta}^{\priv}}
\DeclareMathOperator*{\E}{\mathbb{E}}

\DeclareMathOperator*{\plim}{plim}
\let\vec\relax
\DeclareMathOperator*{\vec}{vec}
\newcommand{\mathbb{R}}{\mathbb{R}}
\newcommand{\mathcal{F}}{\mathcal{F}}
\renewcommand{\P}{\mathbb{P}}
\newcommand{\mathbb{Q}}{\mathbb{Q}}
\newcommand{\mathbf{y}}{\mathbf{y}}
\newcommand{\mathbf{x}}{\mathbf{x}}
\renewcommand{\u}{\mathbf{u}}
\newcommand{w}{w}
\newcommand{\rightsquigarrow}{\rightsquigarrow}
\newcommand{\stackrel{\scriptscriptstyle D}{\approx}}{\stackrel{\scriptscriptstyle D}{\approx}}
\newcommand{\stackrel{P}{\to}}{\stackrel{P}{\to}}
\newcommand{\mathcal{N}}{\mathcal{N}}
\newcommand{\mathbf{0}}{\mathbf{0}}
\DeclareMathOperator{\sign}{sign}
\newcommand{\hat{\beta}}{\hat{\beta}}
\newcommand{\widetilde}{\widetilde}
\newcommand{\mathbf{z}}{\mathbf{z}}
\newcommand{\mathcal{A}}{\mathcal{A}}
\newcommand{\mathcal{B}}{\mathcal{B}}

\DeclareMathOperator{\width}{width}

\definecolor{darkspringgreen}{rgb}{0.09, 0.45, 0.27}
\newcommand{\red}[1]{\textcolor{red}{ #1 }}
\newcommand{\blue}[1]{\textcolor{blue}{ #1 }}
\newcommand{\ds}[1]{\textcolor{red}{\small [#1 \bf--DS]}}
\newcommand{\cf}[1]{\textcolor{blue}{\small [#1 \bf--CF]}}
\newcommand{\sfcomment}[1]{\textcolor{brown}{\small [#1 \bf--SF]}}

\newcommand{\edit}[1]{\textcolor{red}{#1}}

\usepackage{array}
\usepackage{multicol}

\renewcommand{\ds}[1]{}
\renewcommand{\cf}[1]{}

\usepackage{ifthen}
\newcommand{1}{1}

\newcommand{\swap}[3][-]{#3#1#2}

\begin{document}

\twocolumn[

\aistatstitle{Parametric Bootstrap for Differentially Private Confidence Intervals}

\aistatsauthor{
  Cecilia Ferrando, \; Shufan Wang, \; Daniel Sheldon\\
  College of Information and Computer Sciences\\
  University of Massachusetts Amherst\\
  \texttt{\{cferrando\},\{shufanwang\},\{sheldon\}@cs.umass.edu} \\
}

\aistatsaddress{ } ]

\author{%
  Cecilia Ferrando, \; Shufan Wang, \; Daniel Sheldon\\
  College of Information and Computer Sciences\\
  University of Massachusetts Amherst\\
  \texttt{\{cferrando\},\{shufanwang\},\{sheldon\}@cs.umass.edu} \\
}

\begin{abstract}
The goal of this paper is to develop a practical and general-purpose approach to construct confidence intervals for differentially private parametric estimation. We find that the parametric bootstrap is a simple and effective solution. It cleanly reasons about variability of both the data sample and the randomized privacy mechanism and applies ``out of the box'' to a wide class of private estimation routines. It can also help correct bias caused by clipping data to limit sensitivity. We prove that the parametric bootstrap gives consistent confidence intervals in two broadly relevant settings, including a novel adaptation to linear regression that avoids accessing the covariate data multiple times. We demonstrate its effectiveness for a variety of estimators, and find that it provides confidence intervals with good coverage even at modest sample sizes and performs better than alternative approaches.
\end{abstract}

\section{Introduction}\label{introduction}

Differential privacy provides a rubric for drawing inferences from data sets without compromising the privacy of individuals.

This paper is about privately constructing confidence intervals. In the non-private case, approximate methods based on asymptotic normality or the bootstrap~\citep{efron1979bootstrap} apply to a wide range of models and are very widely used in practice.
In the private case, such ``swiss army knife'' methods are hard to find. The situation is complicated by the fact that private estimation procedures are necessarily randomized, which leads to a distinct source of randomness (``privacy noise'') in addition to the usual random draw of a finite sample from a population (``sampling noise'').
We find experimentally that asymptotic methods are significantly less effective in private settings, due to privacy noise that becomes negligible only for for very large sample sizes (Section~\ref{sec:experiments}). Bootstrap approaches face the challenge of incurring privacy costs by accessing the data many times~\citep{brawner2018bootstrap}.

This paper advocates using the parametric bootstrap as a simple and effective method to construct confidence intervals for private statistical estimation. The parametric bootstrap resamples data sets from an estimated parametric model to approximate the distribution of the estimator. It is algorithmically simple, can be used with essentially any private estimator, and cleanly reasons about both sampling noise and privacy noise. Unlike the traditional bootstrap, it is based on post-processing and avoids accessing the data many times, so it often has little or no privacy burden. By reasoning about the distribution of a finite sample, it makes fewer assumptions than purely asymptotic methods and significantly mitigates the problem of non-negligible privacy noise. The parametric bootstrap can also help correct bias in private estimation caused by artificially bounding data to limit sensitivity.

We first introduce the parametric bootstrap and discuss its application to private estimation, including methods to construct confidence intervals and correct bias. We then review parametric bootstrap theory, and apply the parametric bootstrap to obtain provably consistent confidence intervals in two private estimation settings---exponential families and linear regression sufficient statistic perturbation (SSP)---as well as an empirical demonstration for the ``one posterior sample'' (OPS) method~\citep{Wang:2015aa,foulds2016theory,ZZhang:2016aa}. These demonstrate the broad applicability the parametric bootstrap to private estimation.

One limitation of the parametric bootstrap is the restriction to fully parametric estimation. For example, it doesn't apply directly to regression problems that do not have a parametric model for covariates, and may not be appropriate very complex data.
In our linear regression application, we contribute a novel hybrid bootstrap approach to circumvent this limitation; the resulting method is easy to use and simultaneously estimates regression coefficients and constructs confidence intervals with good coverage properties. A second limitation is computational cost, which scales with the data size. For small or medium data sets, the cost is likely manageable. For very large ones, cheap asymptotic methods will often be adequate (see Section~\ref{sec:experiments}; for exponential families and linear regression with sufficient statistic perturbation, the asymptotic distributions are a relatively simple byproduct of our bootstrap theory). However, it is unknown in general how large data must be for asymptotic methods to perform well.

\section{Background}\label{background}

Differential privacy is a formal definition to capture the notion that, to maintain privacy, the output of an algorithm should remain nearly unchanged if the data of one individual changes. Say that two data sets $X$ and $X'$ of size $n$ are \emph{neighbors} if they differ in exactly one data record.

\begin{definition}[Differential privacy, \citealt{dwork2006calibrating}]
A randomized algorithm $\mathcal{A}$ satisfies $\epsilon$-differential
privacy ($\epsilon$-DP) if, for neighboring data sets $X$ and $X'$, and any subset $O \subseteq \text{Range}(\mathcal{A})$,
$$ \Pr[\mathcal{A}(X) \in O] \leq \exp(\epsilon) \Pr[\mathcal{A}(X') \in O]. $$
\end{definition}

One common way to achieve differential privacy is by injecting calibrated noise onto the statistics computed from the data. Let $f$ be any function that maps data sets to $\mathbb{R}^{d}$. The magnitude of noise required to privatize the computation of $f$ depends on its \emph{sensitivity}.

\begin{definition}[Sensitivity, \citealt{dwork2006calibrating}]
The sensitivity of a function f is $$\Delta f = \max_{X,X^{'}} \Vert  f(X) - f(X^{'}) \Vert_{1}$$ where $X,X^{'}$ are any two neighboring data sets.
\end{definition}

When $f$ is additive, it is straightforward to bound its sensitivity (proof in Appendix~\ref{sec:proofs}):
\begin{restatable}{claim}{additive}
\label{claim:additive}
Suppose $X = (x_1, \ldots, x_n)$ and $f(X) = \sum_{i=1}^n g(x_i)$ where $g$ maps data points to $\mathbb{R}^m$. Let $\width (g_j) = \max_{x} g_j(x) - \min_x g_j(x)$ where $x$ ranges over the data domain. Then $\Delta f \leq \sum_{j=1}^m \width(g_j)$, which is a constant independent of $n$.
\end{restatable}

Many algorithms satisfy differential privacy by using the Laplace mechanism.
\begin{definition}[Laplace mechanism, \citealt{dwork2006calibrating}]
Given a function $f$ that maps data sets to $\mathbb{R}^{m}$, the Laplace mechanism outputs the random variable $\mathcal{L}(X) \sim \Lap(f(X), \Delta f / \epsilon)$
from the Laplace
distribution, which has density $\Lap(z; u, b) = (2b)^{-m} \exp ( - \left\Vert  z - u \right\Vert_{1} / b)$. This corresponds to
adding zero-mean independent noise $u_i \sim \Lap(0, \Delta f / \epsilon)$ to each component of $f(X)$.
\end{definition}

\section{Parametric Bootstrap}\label{sec:bootstrap}
\begin{wrapfigure}{R}{0.32\textwidth}
\vspace{-23pt}
\begin{minipage}{0.32\textwidth}
\begin{algorithm}[H]

  \caption{Parametric Bootstrap\label{alg:bootstrap}}
  \begin{algorithmic}[1]

  \State \textbf{Input:} $x_{1:n}$, $B$, estimator $\mathcal{A}$
  \State $\hat{\theta}, \hat \tau \leftarrow \mathcal{A}(x_{1:n})$
    \For {$b$ from $1$ to $B$}
        \State $x^*_1, \ldots, x^*_n \sim P_{\hat{\theta}}$
        \State $\hat\theta^{*b}, \hat\tau^{*b} \leftarrow \mathcal{A}(x^*_{1:n})$
    \EndFor
    \State \Return{$\hat{\tau}, \big(\hat\tau^{*1}, \ldots, \hat\tau^{*B}\big)$}
\end{algorithmic}
\end{algorithm}
\end{minipage}
\vspace{-8pt}
\end{wrapfigure}
We consider the standard setup of parametric statistical inference, where a data sample $x_{1:n} = (x_1, \ldots, x_n)$ is observed and each $x_i$ is assumed to be drawn independently from a distribution $P_\theta$ in the family $\{P_\theta: \theta \in \Theta\}$ with unknown $\theta$.

The goal is to estimate some population parameter $\tau = \tau(\theta)$, the \emph{estimation target}, via an estimator $\hat\tau = \hat\tau(x_{1:n})$.\footnote{We use a hat on variables that are functions of the data and therefore random.} We also seek a $1-\alpha$ confidence interval for $\tau$, that is, an interval $[\hat a_n, \hat b_n]$ such that $\P_\theta( \hat a_n \leq \tau \leq \hat b_n) \approx 1-\alpha$, where $\P_\theta$ is the probability measure over the full sample when the true parameter is~$\theta$.
We will require $\hat\tau$ and $[\hat a, \hat b]$ to be differentially private. Our primary focus is not designing private estimators $\hat\tau$, but designing methods to construct private confidence intervals $[\hat a, \hat b]$ that can be used for many estimators and have little additional privacy burden.

The parametric bootstrap is a simple way to approximate the distribution of $\hat\tau$ for confidence intervals and other purposes. It is a variant of Efron's bootstrap~\citep{efron1979bootstrap,efron1981anonparametric, efron1981bnonparametric, efron1986bootstrap}, which runs an estimator many times on simulated data sets whose distribution approximates the original data. In the parametric bootstrap, data sets are simulated from $P_{\hat\theta}$, the parametric distribution with estimated parameter $\hat\theta$.\footnote{In the non-parametric bootstrap, data sets are simulated from the empirical distribution of $x_{1:n}$.} The procedure is shown in Algorithm~\ref{alg:bootstrap}, where $\mathcal{A}$ is an algorithm that computes the estimates $\hat\theta$ and $\hat\tau$ from the data. A simple case is when $\hat\tau(x_{1:n}) = \tau\big(\hat\theta(x_{1:n})\big)$ but in general these may be estimated separately.

The parametric bootstrap is highly compatible with differential privacy. The data is only accessed in Line 2, so the only requirement is that $\mathcal{A}$ be differentially private (which necessitates it is randomized). The remaining steps are post-processing and incur no additional privacy cost. The simulation cleanly handles reasoning about both data variability (Line 4) and randomness in the estimator (Line 5). When the estimation target is $\theta$, we have $\hat\tau = \hat\theta$, and the procedure incurs no privacy cost beyond that of the private estimator $\hat\tau$. In other cases, additional privacy budget is required to estimate the full vector $\theta$ including nuisance parameters; an example is estimating the mean of a Gaussian with unknown variance~\citep{du2020differentially}.

\subsection{Confidence Intervals and Bias Correction}\label{CIsBiasCorrection}

\begin{figure}[h]
  \centering
\begin{minipage}{0.45\textwidth}
\vspace{-20pt}
\hspace{-15pt}
\begin{table}[H]

  \scriptsize
  \setlength{\tabcolsep}{2.8pt}
  \renewcommand{\arraystretch}{1.2}
  \newcolumntype{L}{>{\raggedright\arraybackslash}m{3cm}}
  \caption{ Bootstrap confidence intervals. $\hat{\xi}_{\gamma}$ is the $1-\gamma$ quantile of $\hat\tau^* - \hat\tau$ or $(\hat\tau^* - \hat\tau)/\hat\sigma^*$, and $\hat{\zeta}_\gamma$ is the $1-\gamma$ quantile of $\hat\tau^*$. \label{table:CI-methods}}

\begin{tabular}{@{}lccc@{}}
  \toprule
  Interval       & Target & Target interval & $\tau$ interval \\
  \midrule
  Pivotal  &
  $\hat\tau - \tau$ &
  $[\hat\xi_{1-\frac{\alpha}{2}}, \hat\xi_{\frac{\alpha}{2}}]$ &
  $[\hat\tau - \hat\xi_{\frac{\alpha}{2}},\ \hat\tau - \hat\xi_{1-\frac{\alpha}{2}}]$ \\
  %
  \midrule
  Studentized\\pivotal &
  $\frac{(\hat\tau - \tau)}{\hat\sigma}$ &
  $[\hat\xi_{1-\frac{\alpha}{2}}, \hat\xi_{\frac{\alpha}{2}}]$ &
  $[\hat\tau - \hat\xi_{\frac{\alpha}{2}} \hat\sigma,\ \hat\tau - \hat\xi_{1-\frac{\alpha}{2}}\hat\sigma]$
  \\
  %
  \midrule
  Efron's\\percentile &
  $\hat\tau$ &
  $[\hat\zeta_{1-\frac{\alpha}{2}},\ \hat\zeta_{\frac{\alpha}{2}}]$ &
  $[\hat\zeta_{1-\frac{\alpha}{2}},\ \hat\zeta_{\frac{\alpha}{2}}]$ \\
  %

  \bottomrule
\end{tabular}
\end{table}
\end{minipage}
\end{figure}

There are several well known methods to compute confidence intervals from bootstrap replicates. Three are listed in Table~\ref{table:CI-methods}; note that names are inconsistent in the literature.\footnote{Our mathematical presentation follows~\citet{van2000asymptotic}, but names  follow~\citet{wasserman2006all}. The names ``pivotal'' and ``studentized pivotal'' and are descriptive and avoid the confusion of ``percentile interval'' sometimes referring to the pivotal interval and other times to Efron's percentile interval. The possessive "Efron's"~\citep{van2000asymptotic} clarifies that we use Efron's definition of ``percentile interval''~\citep[e.g.,][]{efron2016computer}.}
The general principle is to treat the pair $(\hat\tau^*, \hat\tau)$ analogously to $(\hat\tau, \tau)$ to approximate the distribution of the latter. The intervals differ according to what function of $(\hat\tau, \tau)$ they target. A simple example is to approximate the ``pivot'' $\hat\tau - \tau$ by $\hat\tau^* - \hat\tau$, which leads to the pivotal interval. To construct it, we estimate the $1-\gamma$ quantile of $\hat\tau^* - \hat\tau$ as the $1-\gamma$ quantile of the bootstrap replicates $(\hat\tau^{*1}-\hat\tau, \ldots, \hat\tau^{*B}-\hat\tau)$. The number of replicates controls the error introduced by this step. This error is usually ignored theoretically because, in principle, it can be reduced arbitrarily with enough computation, and it can be controlled well in practice. The studentized pivotal interval targets $(\hat\tau - \tau)/\hat\sigma$ instead, where $\hat\sigma$ is a standard error estimate of the main procedure; it can converge faster than the pivotal interval~\citep{wasserman2006all}. Efron's percentile interval targets $\hat\tau$ directly and, while simple, its logic is less obvious; it can also be viewed as targeting the pivot $\hat\tau - \tau$ with a ``reversed'' interval, which is how theoretical properties are shown. By approximating $\hat\tau - \tau$ by $\hat\tau^* - \hat\tau$, we can also estimate the bias of $\hat\tau$. This leads to a simple bias corrected estimator $\hat\tau_{\text{bc}}$:
$$
\widehat{\text{bias}} = \mathbb{E}[\hat\tau^* - \hat\tau], \quad \hat{\tau}_{\text{bc}} \leftarrow \hat \tau - \widehat{\text{bias}}.
$$
Similar to the quantiles above, $\mathbb{E}[\hat\tau^* - \hat\tau]$ is estimated as the sample mean over bootstrap replicates.

\subsection{Significance and Connection to Other Resampling Methods for Private Estimation}\label{significance}

The parametric bootstrap can be applied to any parametric estimation problem, a wide range of private estimators, is very accurate in practice, and has little or no additional cost in terms of privacy budget or algorithm development.
These make it an excellent default choice (to our knowledge, the best) for constructing private confidence intervals for any parametric estimation problem with small to medium data sets.

That such a simple and effective choice is available is not articulated in the literature. Two prior works use methods that can be viewed as the parametric bootstrap, but do not discuss the classical procedure and its wide ranging applications, or different facets of bootstrap methodology such as techniques for constructing confidence intervals and bootstrap theory. Specifically, the simulation approach of \citet{du2020differentially} for Gaussian mean confidence intervals is equivalent to the parametric bootstrap with a non-standard variant of Efron's percentile intervals, and performed very well empirically. In their application to independence testing, \citet{gaboardi2016differentially} approximate the distribution of a private test statistic by simulating data from a null model after privately estimating its parameters; this can be viewed as an application of the parametric bootstrap to the null model.

Several other works use resampling techniques that resemble the parametric bootstrap for a similar, but conceptually distinct, purpose~\citep{d2015differential,wang2019differentially,evans2019statistically}. A typical setup is when $\hat\tau = \tau' + \eta$, where $\tau'$ is a non-private estimator and $\eta$ is noise added for privacy. Standard asymptotics are used to approximate $\sqrt{n}(\tau' - \tau)$ as $\mathcal{N}(0, \hat\sigma)$, where $\hat\sigma$ is a (private) standard error estimate for $\tau'$. For the private estimator, this gives $\sqrt{n}(\hat\tau-\tau) \approx \mathcal{N}(0, \hat\sigma) + \sqrt{n}\eta$. Because $\eta$ has known distribution, Monte Carlo sampling can be used to draw samples from $\mathcal{N}(0, \hat\sigma) + \sqrt{n}\eta$ for computing confidence intervals or standard errors. The key distinction is that \emph{standard} asymptotics are used to approximate the distribution of $\tau'$, which captures all variability due to the data, and sampling is used only to combine that distribution with the privacy noise distribution. In contrast, the key feature of a bootstrap method is that it resamples data sets to reason about estimator variability due the random data, and thereby avoids standard asymptotics. This technique also does not apply when the privacy mechanism is more complicated than adding additive noise to a non-private estimate (cf. the OPS example of Sec.~\ref{sec:OPS}).

\section{Bootstrap Theory}\label{theory}

This section gives general results we can use to argue correctness of bootstrap confidence intervals in private settings. We give a general notion of ``bootstrap'' estimator that covers different  resampling methods. Let $(\Omega, \mathcal{F}, \P_\theta)$ be the probability space for $x_1, x_2, \ldots \sim P_\theta$ and $\eta_1, \eta_2, \ldots$ where, for a given $n$, the data is $x_{1:n}$ and $\eta_n$ captures any other randomness used in the privacy mechanism or estimator; we refer to this as the ``outer'' probability space. A bootstrap estimator is defined in terms of a random experiment over an ``inner'' probability space conditional on $\omega \in \Omega$ and $n$. Let $\P_n^*(\cdot | \omega)$ be a Markov kernel defining this space. The traditional bootstrap uses $\P_n^*(\cdot | \omega) = \hat{P}^n$ with $\hat{P}(dx) = \frac{1}{n} \sum_{i} \delta_{x_i}(dx)$; the parametric bootstrap uses $\hat{P} = P_{\hat \theta}$ instead. Our hybrid model in Section~\ref{linearregression} uses a custom resampling method, which gives a custom measure $\P_n^*(\cdot | \omega)$.

For our purposes, a bootstrap estimator of a parameter $\tau(\theta)$ is a random variable $\hat\tau^*_n$ in the inner probability space that simulates the parameter estimate $\hat\tau_n$ of the ``main'' procedure. Typically, the bootstrap estimator arises from running the main procedure on resampled data. That is, if $\hat\tau_n = T_n(\omega)$, then $\hat\tau^* = T_n(\omega^*)$ with $\omega^*  \sim \P_n^*(\cdot \mid \omega)$. Our hybrid OLS bootstrap will deviate slightly from this pattern.

\subsection{Consistency}
Bootstrap ``success'' has to do with the asymptotic distributions of the (approximate) pivot $\sqrt{n}(\hat\tau_n - \tau)$ and its bootstrapped counterpart $\sqrt{n}(\hat\tau^*_n - \hat\tau_n)$. For studentized intervals, the pivot $(\hat\tau_n - \tau)/\hat \sigma_n$ is used instead, where $\hat\sigma_n$ is a standard error estimate of the main procedure; theory for this case is a straightforward extension if $\hat\sigma_n \to \sigma(\theta)$ in $\P_\theta$-probability~\citep{van2000asymptotic,beran1997diagnosing}.

\begin{definition}
  \label{def:consistency}
  The bootstrap estimator $\hat\tau^*_n$ is \emph{consistent} if
\begin{equation}
\begin{split}
\sup_x
\Big|
\P_n^*\Big(
   \sqrt{n}(\hat\tau_n^* - \hat\tau_n) \leq t \mid \omega
\Big)
- \\-
\P_\theta\Big(
   \sqrt{n}(\hat\tau_n - \hat\tau) \leq t
   \Big) \Big|
\stackrel{P}{\to} 0
\end{split}
\end{equation}
with convergence in $\P_\theta$-probability.
\end{definition}
This says that the Kolmogorov-Smirnov distance between the distribution of the pivot and the conditional distribution of the bootstrapped pivot converges to zero, in probability over $\omega$.

In most cases $\sqrt{n}(\hat \tau_n - \tau) \rightsquigarrow T$ for a continuous random variable $T$. In this case it is enough for the bootstrapped pivot to converge to the correct limit distribution.
\begin{lemma}[\citealt{van2000asymptotic}, Eq. (23.2)]
  \label{lem:consistency}
  Suppose $\sqrt{n}(\hat \tau_n - \tau) \rightsquigarrow T$ for a random variable $T$ with continuous distribution function $F$. Then, $\hat\tau^*_n$ is consistent if and only if, for all $t$,
  $$
  \P^*_n \Big(\sqrt{n}(\hat\tau^*_n - \hat\tau_n) \leq t \mid \omega \Big) \overset{P}{\to} F(t).
  $$
\end{lemma}

Consistency is also preserved under continuous mappings: if $\hat\tau^*_n$ is consistent relative to $\sqrt{n}(\hat\tau_n - \tau) \rightsquigarrow T$ and $g$ is continuous, then $g(\hat\tau^*_n)$ is consistent relative to $\sqrt{n}\big(g(\hat\tau_n) - g(\tau)\big)$~\citep{beran1997diagnosing}. In our applications we will show consistency of a bootstrap estimator $\hat\theta^*_n$ for the full parameter vector $\theta$, which implies consistency for continuous functions of $\theta$; a simple application is selecting one entry and constructing a confidence interval.

\subsection{Confidence interval consistency}

Bootstrap consistency implies consistent confidence intervals. The confidence interval $[\hat{a}_n, \hat b_n]$ for $\tau=\tau(\theta)$ is (conservatively) asymptotically consistent at level $1 - \alpha$ if, for all $\theta$,
\begin{equation}
  \liminf_{n \to \infty} \P_\theta \left(\hat a_{n} \leq \tau \leq \hat b_{n}\right) \geq 1 - \alpha.
\end{equation}

\begin{lemma}[\citealt{van2000asymptotic}, Lemma 23.3]
Suppose $\sqrt{n}(\tau_n - \tau) \rightsquigarrow T$ for a random variable $T$ with continuous distribution function and $\tau^*_n$ is consistent. Then the pivotal intervals are consistent, and, if $T$ is symmetrically distributed around zero, then Efron's percentile intervals are consistent. When the analogous conditions hold for the studentized pivot $(\hat\tau_n - \tau)/\hat\sigma_n$, studentized intervals are consistent.
\end{lemma}

\subsection{Parametric bootstrap consistency}

\citet{beran1997diagnosing} showed that asymptotic equivariance of the main estimator guarantees consistency of the parametric bootstrap. Let $H_n(\theta)$ be the distribution of $\sqrt{n}(\hat\tau_n - \tau(\theta))$ under $\P_\theta$.
\begin{definition}[Asymptotic equivariance,~\citealt{beran1997diagnosing}]\label{definition:AE}
The estimator $\hat\tau_n$ is asymptotically equivariant if $H_n(\theta + h_n/\sqrt{n})$ converges to a limiting distribution $H(\theta)$ for all convergent sequences $h_n$ and all~$\theta$.
\end{definition}
\begin{restatable}[Parametric bootstrap consistency]{theorem}{PBconsistency}\label{PBconsistency}
Suppose $\sqrt{n}(\hat\theta_n - \theta) \rightsquigarrow J(\theta)$ and $\hat\tau_n$ is asympotitcally equivariant with continuous limiting distribution $H(\theta)$. Then the parametric bootstrap estimator $\hat\tau^*_n$ is consistent.
\end{restatable}
All proofs are provided in the appendix.
Furthermore, under reasonably general conditions, the reverse implication is true, with bootstrap failures occurring precisely at those parameter values $\theta_0$ for which asymptotic equivariance does not hold~\citep{beran1997diagnosing}.

\section{Applications}

We apply the parametric bootstrap to three private estimation settings: (i)~exponential families with sufficient statistic perturbation (SSP), (ii)~linear regression with SSP, (iii)~the ``one posterior sample'' (OPS) estimator.

\paragraph{Exponential Families}\label{exponentialfamilies}

A family of distributions is an \emph{exponential family} if $P_\theta$ has a density of the form:
$$
\vspace{-2pt}
p(x; \theta) = h(x) \exp(\theta^T T(x) - A(\theta))
$$
where $h(x)$ is a base measure, $\theta$ is the natural parameter, $T$ is the sufficient statistic function, and $A$ is the log-partition function. Define the log-likelihood function of an exponential family as $\ell(\theta; x) = \log p(x; \theta) - \log h(x) = \theta^T T(x) - A(\theta)$.
The constant term $\log h(x)$ does not affect parameter estimation and is subtracted for convenience.
For a sample $x_{1:n}$, let $T(x_{1:n}) = \sum_{i=1}^n T(x_i)$. The log-likelihood of the sample is
$$
\ell(\theta; x_{1:n}) = \theta^T T(x_{1:n}) - n A(\theta) := f\big(\theta; T(x_{1:n}) \big),
$$
which depends on the data only through the sufficient statistic $T(x_{1:n})$. The maximum-likelihood estimator (MLE) is
$\hat{\theta} = \argmax_\theta f\big(\theta; T(x_{1:n})\big)$.

A simple way to create a private estimator is sufficient statistic perturbation (SSP); that is, to privatize the sufficient statistics using an elementary privacy mechanism such as the Laplace or Gaussian mechanism prior to solving the MLE problem.
SSP is a natural choice because $T(x_{1:n})$ is a compact summary of the data and has sensitivity that is easy to analyze, and it often works well in practice~\citep{bernstein2018differentially,foulds2016theory}.
Specifically, it means solving
\begin{equation}
\hat{\theta} = \argmax_\theta f\big(\theta, T(x_{1:n}) + w\big) \tag{SSP-MLE}
\label{SSPMLE}
\end{equation}
where $w$ is a suitable noise vector. This problem has closed form solutions for many exponential families and standard numerical routines apply to others.
For the Laplace mechanism, $w_j \sim \Lap(\frac{\Delta}{\epsilon})$
\ds{Dropped the factor of 2. BTW, what do we do in our code?}\cf{we use the range of a surrogate sample of size 1000 (the 2x was an old strategy to stay conservative, we dropped in the code long ago)}
for all $j$, where $\Delta = \sum_{j} \width(T_j)$ is an upper bound on the $L_1$ sensitivity of $T(x_{1:n})$ by Claim~\ref{claim:additive}. If $\width(T_j)$ is not known or is unbounded, the analyst must supply bounds and guarantee they are met, e.g., by discarding data points that don't meet the bounds, or clamping them to the bounded interval.

\begin{restatable}{theorem}{MLEconsistency}\label{MLEconsistency}
Let $\hat{\theta}_n$ be the solution to the (SSP-MLE) optimization problem for a sample $x_{1:n}$ from an exponential family model that satisfies the regularity conditions given in \citet[Section 4.4.2]{davison2003statistical}. Then $\sqrt{n}(\hat{\theta}_n - \theta)$ is asymptotically equivariant with limiting distribution $\mathcal{N}(0, I(\theta)^{-1})$, where $I(\theta) = \nabla^2 A(\theta)$ is the Fisher information. This implies consistency of the parametric bootstrap estimator $\hat\theta^*_n$.
\end{restatable}
\ds{Give reference to specific appendix?}

\paragraph{Linear Regression}\label{linearregression}
\label{sec:OLS}

We consider a linear regression model where we are given $n$ pairs\footnote{We use boldface for vectors as needed to distinguish from scalar quantities.}  $(\mathbf{x}_i, y_i)$ with $\mathbf{x}_i \in \mathbb{R}^p$ and $y_i \in \mathbb{R}$ assumed to be generated as
$
y_i = \beta^T \mathbf{x}_i + u_i,
$
where the errors $u_i$ are i.i.d., independent of $\mathbf{x}_i$, zero-mean, and have finite variance $\sigma^2$, and the $\mathbf{x}_i$ are i.i.d. with $\E[\mathbf{x} \mathbf{x}^T] = Q$. We wish to estimate the regression coefficients $\beta \in \mathbb{R}^p$. Let $X \in \mathbb{R}^{n \times p}$ be the matrix with $i$th row equal to $\mathbf{x}_i^T$ and $\mathbf{y}, \u \in \mathbb{R}^N$ be the vectors with $i$th entries $y_i$ and $u_i$, respectively. The \emph{ordinary least squares} (OLS) estimator is:
\begin{equation}
  \hat \beta = (X^{T}X)^{-1} X^{T} \mathbf{y}.
  \label{eq:OLS}
\end{equation}

Like the MLE in exponential families, Eq.~\eqref{eq:OLS} depends on the data only through sufficient statistics $X^T X$ and $X^T \mathbf{y}$, and
SSP is a simple way to privatize the estimator that works very well in practice~\citep{wang2018revisiting}. The privatized estimator is
\begin{equation} \label{eq:privnorm}
    \hat{\beta} = (X^{T}X + V)^{-1} (X^{T} \mathbf{y} + w),
    \tag{SSP-OLS}
\end{equation}
where $V \in \mathbb{R}^{p \times p}$ and $w \in \mathbb{R}^p$ are additive noise variables drawn from distributions $P_V$ and $P_w$ to ensure privacy.

For the Laplace mechanism, we use \ds{This paragraph could go to appendix, maybe along with details of $\hat\sigma^2$? We do need to mention $\hat\sigma^2$; the properties we need are that it is private and consistent.}
\begin{align}
V_{jk} &\sim \Lap(0, \Delta_{V} /\epsilon_1) \text{ for } j \leq k, \text{ and } V_{kj} = V_{jk}, \label{Vterm} \\
w_j &\sim \Lap(0, \Delta_{w}/ \epsilon_2) \text{ for all $j$}, \label{Wterm}
\end{align}
where $\Delta_{V}$ and $\Delta_{w}$ bound the $L_1$ sensitivity of $V$ and $w$, respectively. The result is $(\epsilon_1 + \epsilon_2)$-DP. Because $X^TX = \sum_{i=1}^n \mathbf{x}_i \mathbf{x}_i^T$ and $X^T \mathbf{y} = \sum_{i=1}^n \mathbf{x}_i y_i$ are additive, we can take $\Delta_V = \sum_{j \leq k} \width(x_j)\cdot\width(x_k)$ and $\Delta_w = \sum_j \width(x_j)\cdot\width(y)$, where $\width(x_j)$ and $\width(y)$ are widths of the $j$th feature and response variable, respectively, which are enforced by the modeler.

For confidence intervals, we will also need a private estimate of $\sigma^2$: let $\hat{\sigma}^2 = (n-p)^{-1}\sum_{i=1}^n (y_i - \hat{\beta}^T \mathbf{x}_i)^2 + \Lap(0, \Delta_z / \epsilon_3)$ where $\Delta_z = \width((y - \hat{\beta}^T \mathbf{x})^2)$. \ds{TODO} The released values for SSP-OLS are then $(X^TX + V, X^T \mathbf{y} + w, \hat\beta, \hat\sigma^2)$, which satisfy $(\epsilon_1 + \epsilon_2 + \epsilon_3)$-DP.

\paragraph{Limitations of parametric bootstrap for private regression}
The parametric bootstrap is more difficult to apply to regression problems in a private setting due to the covariates. It is typical to bootstrap conditioned on $X$, which means simulating new response variables $\mathbf{y}$ from a parametric distribution $p(\mathbf{y} | X; \hat{\beta}, \hat{\sigma}^2)$, where $\hat\beta$ and $\hat\sigma^2$ are (privately) estimated parameters, and a fully parametric distribution $p(u; \sigma^2)$ is assumed for errors. A bootstrap replicate would look like $\hat\beta^* = (X^TX)^{-1}X^T\mathbf{y}^*$ with $\mathbf{y}^* = X\hat\beta + \u^*$ and $\u^*$ simulated form the error distribution.
The challenge is that $X$ is accessed to generate each replicate, so to make it differentially private would require additional randomization and consume privacy budget. An alternative would be to posit a model $p(\mathbf{x}; \theta)$ and perform the parametric bootstrap with respect to a joint model $p(\mathbf{x}, y ; \theta, \beta, \sigma^2)$, but the additional demand to model covariates is unappealing in a regression context.

\paragraph{Hybrid parametric bootstrap for OLS}
We propose a novel hybrid approach that avoids the need to repeatedly access covariate data or to model the covariate or error distributions explicitly. Conceptually, we use the part of the standard asymptotic analysis that ``works well'' to approximate the relevant statistics of the covariate data, and use the parametric bootstrap to deal with the noise added for privacy.
Following standard analysis for OLS, we can substitute $\mathbf{y} = X \beta + \u$ in \eqref{eq:privnorm} and scale terms to get:
\begin{equation}
\begin{split}
\label{eq:beta_hat}
\sqrt{n}\hat{\beta}_n  &=
\sqrt{n}\left(\check Q_n + \frac{1}{n}V \right)^{-1} \check Q_n  \beta  +
\\ &+ \left(\check Q_n + \frac{1}{n}V \right)^{-1} \left(\blue{\frac{1}{\sqrt{n}} X^T \mathbf{u}} + \frac{1}{\sqrt{n}} w \right), \\ \,\, \check Q_n &= \frac{1}{n} X^T X.
\end{split}
\end{equation}

This expression is instructive to see the different sources of randomness that contribute to the variability of $\hat\beta_n$:
the terms $\check Q_n = \frac{1}{n}X^TX$ and $\frac{1}{\sqrt{n}}X^T \u$ are due to data variability, and $\frac{1}{n}V$ and $\frac{1}{\sqrt{n}}w$ are due to privacy.
We form a bootstrap estimator $\hat\beta^*_n$ that treats $(\hat\beta_n, \beta)$ analogously to $(\hat\beta^*_n, \hat\beta_n)$ and simulates the different sources of variability using the best available information about their distributions:
\begin{equation}
\begin{split}
\label{eq:beta_hat_star}
  \sqrt{n}\hat{\beta}^*_n &= \sqrt{n}\left(\hat{Q}_n + \frac{1}{n}V^* \right)^{-1} \hat Q_n \hat{\beta}_n +
  \\ &+ \left(\hat{Q}_n + \frac{1}{n}V^* \right)^{-1} \left(\blue{Z^*_n} + \frac{1}{\sqrt{n}} w^* \right), \\
  \notag \blue{Z_n^*} &\sim \mathcal{N}(0, \hat\sigma^2_n \hat Q_n), \quad V^* \sim P_V, \quad w^* \sim P_w, \\
  \notag \hat Q_n &= \frac{1}{n}X^TX + \frac{1}{n}V.
 \end{split}
\end{equation}

All privacy terms in Eq.~\eqref{eq:beta_hat} are simulated from their exact distributions in Eq.~\eqref{eq:beta_hat_star}. The variables $\check Q_n$ and $\hat Q_n$ represent approximations of $Q = \E[\mathbf{x}\x^T]$ available to the corresponding estimator. Both quantities converge in probability to $Q$. Our choice not to simulate variability in these estimates due to the covariates is analogous to the ``fixed $X$'' bootstrap strategy for regression problems~\citep{fox2002r}; we \emph{do} simulate the variability due to privacy noise added to the estimates. The blue terms represent contributions to estimator variability due to interactions between covariates and unobserved noise variables. In a traditional bootstrap, we might simulate this term in Eq.~\eqref{eq:beta_hat_star} as $\frac{1}{\sqrt{n}}X^T \u^*$ where $\u^*$ are simulated errors, but, as described above, we do not wish to access $X$ within the bootstrap procedure. Instead, because we know $\frac{1}{\sqrt{n}}X^T \u \rightsquigarrow \mathcal{N}(0, \sigma^2 Q)$ by the central limit theorem,\footnote{This is a standard result of OLS asymptotics and is expected to be accurate for modest sample sizes.} and because $\sigma^2$ and $Q$ are estimable, we simulate this term directly from the normal distribution with estimated parameters.

\begin{restatable}{theorem}{OLS}\label{OLS}
  \label{thm:OLS}
  The private estimator satisfies $\sqrt{n}(\hat\beta_n - \beta) \rightsquigarrow \mathcal{N}(0, \sigma^2 Q^{-1})$ and the bootstrap estimator $\hat\beta^*_n$ is consistent in the sense of Lemma~\ref{lem:consistency}.
\end{restatable}

\paragraph{OPS}
\label{sec:OPS}
\cite{Dimitrakakis:2014aa}, \cite{Wang:2015aa} and \cite{foulds2016theory} used the idea of sampling from a Bayesian posterior distribution to obtain a differentially private point estimate.
One Posterior Sampling (OPS), which releases one
sample from the posterior, is a special case of the exponential mechanism, and the corresponding estimator is near-optimal for parametric learning \citep{Wang:2015aa}.
The parametric bootstrap applies easily to OPS estimators and produces well calibrated intervals (Figure~\ref{fig:fig2}).
We expect the asymptotic analysis of \cite{Wang:2015aa} can be adapted to prove asymptotic equivariance, and hence parametric bootstrap consistency, for OPS, but do not give a formal proof.

\section{Related work}
\label{sec:related}
A number of prior works have studied private confidence intervals for different models~\citep{d2015differential,karwa2017finite,Sheffet:2017aa,Barri:2018,gaboardi2018locally,brawner2018bootstrap,du2020differentially}.
\citet{smith2011privacy} showed that a broad class private estimators based on subsample \& aggregate~\citep{nissim2007smooth} are asymptotically normal.
\cite{d2015differential} proposes an algorithm based on subsample \& aggregate to approximate the variance of a private estimator (see Section~\ref{significance}).
The topics of differentially private hypothesis testing \citep{vu2009differential,solea2014differentially,gaboardi2016differentially,couch2019differentially}
and Bayesian inference~\citep{Williams:2010aa,Dimitrakakis:2014aa,Wang:2015aa,foulds2016theory,ZZhang:2016aa,heikkila2017differentially,bernstein2018differentially,bernstein2019differentially} are also related, but the specific considerations differ somewhat from confidence interval construction.
Finding practical and general-purpose algorithms for differentially private confidence intervals has been identified as an important open problem~\citep{opendp}.

The confidence interval approach of \citet{wang2019differentially} applies to any model fit by empirical risk minimization with objective or output perturbation and is similar to the asymptotic methods we compare to in Section~\ref{sec:experiments}.
\citet{evans2019statistically} also give a general-purpose procedure based on subsample \& aggregate (S\&A)~\citep{nissim2007smooth} with normal approximations.
This method also uses S\&A for the point estimates.
We compare to a similar variant of S\&A in Section~\ref{sec:experiments}.
\cite{wang2018statistical} study statistical approximating distributions for differentially private statistics in a general setting.

\citet{brawner2018bootstrap} use the \emph{non-parametric} bootstrap in a privacy context to estimate standard errors ``for free'' (at no additional cost beyond mean estimation) in some settings.
Other methods most similar to the our work on the parametric bootstrap were discussed in more detail in Sec.~\ref{significance}.

Prior methods to construct confidence intervals for private linear regression include \citep{Sheffet:2017aa,Barri:2018}.

\section{Experiments}\label{sec:experiments}
\begin{figure}[h!]
    \centering
    \includegraphics[width=1.02\linewidth]{images/fig1aexpanded.pdf}
    \caption{Observed vs nominal coverage of 95\% CIs for different distributions for different $n$.}
    \label{fig:1a}
\end{figure}

\begin{figure*}
    \centering
    \includegraphics[width=\linewidth]{images/uniquefigurenew.pdf}
    \caption{Left pair. (i) observed vs. nominal coverage for different coverage levels, for a Poisson, n=100, $\epsilon = 0.5$. (ii) same plot, comparing the OPS method \citep{foulds2016theory} and the parametric bootstrap for Bernoulli estimation. Center pair: average CI widths for different $n$ for (i) Gamma and (ii) OLS (other distributions give qualitatively similar results). Width of the private bootstrap CIs approaches that of the public CIs as $n \rightarrow \infty$. Right pair. (i) private and bias-corrected private estimates for a Poisson clamped at varying right-tail thresholds. (ii) same for a Gaussian clamped at $-10$ on the left tail and at varying thresholds on the right tail.}
    \label{fig:fig2}
\end{figure*}

We design synthetic experiments to demonstrate our proposed methods for differentially private confidence interval estimation.

First, we evaluate the performance of private parametric bootstrap CIs vs a baseline method (``Fisher CIs'') based on asymptotic normality of the private estimator and described in more detail below. Performance is measured by how well the coverage of the private CIs matches the nominal coverage. For all models, we also include Fisher CIs of non-private estimators for comparison.

Second, we demonstrate the bias-correction procedure in Sec.~\ref{CIsBiasCorrection} in the case of Gaussian and Poisson distributions with data points clamped to different thresholds, which introduces estimation bias. These results show the effectiveness of the parametric bootstrap at approximating and mitigating the bias of differentially private estimates when sensitivity is bounded by forcing data to take values within given bounds.

Third, we compare parametric bootstrap CIs to another general purpose method to construct confidence intervals based on subsample \& aggregate~\citep{nissim2007smooth,smith2011privacy,d2015differential}.

Finally, the appendix includes additional experiments exploring a broader range of settings and performance metrics. These include: multivariate distributions, the effect of varying $\epsilon$, and measurements of the upper- and lower-tail CI failures. We aim for private CIs to be as tight as possible while providing the correct coverage: in the appendix, we also compare the width of our intervals with that of intervals from existing methods for the specific case of Gaussian mean estimation of known variance.

\newcommand{\myparagraph}[1]{\smallskip\noindent \textbf{#1}\ }

\myparagraph{Baseline: ``Fisher CIs''}
As a byproduct of our consistency analysis we also derive asymptotic normal distributions of the private estimators for both exponential families (Theorem~\ref{MLEconsistency}) and OLS (Theorem~\ref{thm:OLS}).
In each case, we obtain a private, consistent estimate $\hat \sigma^2_j$ of the $j$th diagonal entry of the inverse Fisher information matrix of the private estimator $\hat\theta_n$, and then construct the confidence interval for $\theta_j$ as
\begin{equation}
C_n = \big[\hat{\theta}_{n,j} - z_{\alpha/2} \hat{\sigma}_j, \hat{\theta}_{n,j} + z_{\alpha/2} \hat{\sigma}_j\big],
\label{standardCI}
\end{equation}
where $z_{\gamma}$ is the $1-\gamma$-quantile of the standard normal distribution.
For exponential families, the Fisher information is estimated via plug-in estimation with the private estimator $\hat\theta_n$. For OLS, it is estimated via plugging in private estimates $\hat Q_n = \frac{1}{n} X^TX + \frac{1}{n}V$ and $\hat \sigma^2_n$, which are both released by the SSP mechanism.
For non-private Fisher CIs, we follow similar (and very standard) procedures with non-private estimators.

\myparagraph{Exponential families}
We use synthetic data sets drawn from different exponential family distributions.
Given a family, true parameter $\theta$, and data size $n$, a data set is drawn from $P_{\theta}$. We release private statistics via SSP with the Laplace mechanism. To simulate the modeler's domain knowledge about the data bounds, we draw a \emph{separate} surrogate data set of size 1000 drawn from the same distribution, compute the data range and use it to bound the width of each released statistic. For private estimation, sampled data is clamped to this range. Private $\hat{\theta}$ is computed from the privately released statistics using \ref{SSPMLE}. For the parametric bootstrap CIs, we implement Algorithm \ref{alg:bootstrap} and compute Efron's percentile intervals (see Table~\ref{table:CI-methods}). The output coverage is computed over $T=1000$ trials.

Results are shown in Figures~\ref{fig:1a} and \ref{fig:fig2}.
For the parametric bootstrap, actual coverage closely matches the nominal coverage, even for very small $n$. Coverage of private Fisher CIs is too low until $n$ becomes large, due to the fact that it ignores privacy noise. The bootstrap procedure correctly accounts for the increased uncertainty due to the privacy by enlarging the CIs. The width of the bootstrap intervals approaches the width of the baseline Fisher intervals as $n \rightarrow \infty$. In the appendix, we show that the coverage failures are balanced between left and right tails and examine the effect of increasing $\epsilon$ (which reduces privacy noise and has the same qualitative effect as increasing $n$).

\myparagraph{Linear regression}
We follow a very similar procedure for OLS. Data is generated with $x_j \sim \text{Unif}([-5, 5])$ for all $j$ and errors are $u_i \sim \text{Unif}[-10, 10]$; bounds on $y$ are passed as inputs ($[-150, 150]$) and assumed known. Observed values of $y$ exceeding the given bounds are dropped. These bounds are also used to compute widths for the sensitivity. Private coefficients are estimated with SSP-OLS and bootstrap CIs are constructed via Efron's percentile method. The results are shown in Fig.~\ref{fig:1a}.

\myparagraph{Bias correction}
In the case of distributions with infinite support, one option to bound the sensitivity is to clamp or truncate the data to given bounds. These procedures may induce estimation bias. As discussed in Sec~\ref{CIsBiasCorrection}, the parametric bootstrap can be used to approximate this bias and mitigate it. We demonstrate bias correction on the private estimates and CIs of a Poisson and Gaussian distribution where data is clamped on the right tail at different thresholds (Fig.~\ref{fig:fig2}).

\myparagraph{Comparison with subsample \& aggregate}
We compare the parametric bootstrap CIs with the intervals obtained via a subsample \& aggregate (S\&A) algorithm.
We adapted the S\&A procedure of \cite{d2015differential} for privately estimating standard errors to compute confidence intervals; see Algorithm~\ref{SACI} in the appendix. We compare the accuracy of point estimates and 95\% CIs for the mean of a Gaussian of known variance.  We found that the parametric bootstrap provides more accurate point estimates and better calibrated, tighter CIs than S\&A (Figure~\ref{fig:sa1}).

\begin{figure}[h!]
    \centering
    \includegraphics[scale=0.4]{images/SAPBgaussian.pdf}
    \caption{Point estimates (left), 95\% CI coverage (center) and average CI width (right) of the S\&A method  (see Algorithm~\ref{SACI} in appendix) vs  parametric bootstrap for the mean of a Gaussian of known variance. Settings: $\epsilon = 0.5$, $\theta=0$, $\sigma=1, (x_{min}, x_{max}) = (-20, +20), (L_{min}, L_{max}) = (-10, +10), var_{max} = 50$.}
    \label{fig:sa1}
\end{figure}

\section{Conclusion}
The parametric bootstrap is useful and effective to construct consistent differentially private confidence intervals for broad classes of private estimators, including private linear regression, for which we present a novel adaptation to avoid accessing the covariate data many times. The parametric bootstrap yields confidence intervals with good coverage even at modest sample sizes, and tighter than the ones based on subsample \& aggregate or other general methods. It can be used with any privacy mechanism, and can help mitigate differentially private estimation bias.

\newpage


\section{Introduction}
\label{Sec_Intro}

While diffusion MRI has provided enhanced sensitivity to tissue microstructure \textit{in vivo} by capturing the translational motion of water molecules diffusing in biological tissue,~\citep{LeBihan:1990, LeBihan:1992, Basser:1994, Jones:2010} diffusion-relaxation MRI (drMRI) additionally reports on the local chemical composition of the aqueous phase.~\citep{Zhang:2014,deSantis_relaxometry:2016, deSantis_T1:2016, Hutter:2018,Park_ISMRM:2018,Slator:2019,deAlmeidaMartins_Topgaard:2018,deAlmeidaMartins:2020} For instance, the longitudinal relaxation time $T_1$ informs on molecular dynamics~\citep{Bloembergen:1947} and on the presence of paramagnetic species~\citep{Bloembergen:1948} in simple liquids. \textit{In vivo}, it is mainly determined by cross relaxation, magnetization transfer and spin diffusion with macromolecules in general~\citep{Edzes_Samulski:1977, Halle:2006,Rooney:2007} and myelin lipids in particular,~\citep{Mottershead:2003, Bjarnason:2005, deSantis:2014, Lutti:2014, Stuber:2014} as well as by the interplay between relaxation and diffusion.~\citep{Brownstein_Tarr:1977} However, microstructural studies have been hindered by the fact that the measured drMRI signal is only sensitive to the voxel-averaged diffusion-relaxation profile, with typical cubic-millimeter voxels comprising multiple cell types and the extra-cellular space.~\citep{Stanisz:1997,Norris:2001,Sehy:2002,Minati:2007,Mulkern:2009}

Two strategies were explored to alleviate the lack of specificity of the drMRI signal.
On the one hand, multiple models and signal representations have been developed to relate either the diffusion-$T_2$~\citep{Veraart:2018,Lemberskiy:2018,Ning:2020} or diffusion-$T_1$~\citep{deSantis_relaxometry:2016, deSantis_T1:2016,Andrews_ISMRM:2019} MRI signal to the voxel content. However, models/signal representations rely on compartmental/functional assumptions that may disagree with the underlying tissue microstructure.~\citep{Jelescu:2017,Novikov_on_modeling:2018,Reymbaut_accuracy_precision:2020}
On the other hand, `tensor-valued' diffusion encoding gradient waveforms have enhanced the specificity of the data itself by targeting specific features of the intra-voxel diffusion profile.~\citep{Eriksson:2013,Westin:2014,Eriksson:2015,Westin:2016,Topgaard:2017,Topgaard_dim_rand_walks:2019} Indeed, such measurements typically use four acquisition dimensions - the trace $\mathit{b}$ (size), normalized anisotropy $\mathit{b}_\Delta\in [-0.5,1]$ (shape) and orientation $(\Theta,\Phi)$ of an axisymmetric encoding tensor $\mathbf{b}$~\citep{Basser:1994,Mattiello:1994,Mattiello:1997} - to probe the four dimensions of microscopic axisymmetric diffusion tensors, \textit{i.e.} their isotropic diffusivity $\mathit{D}_\mathrm{iso}$, normalized anisotropy $\mathit{D}_\Delta\in [-0.5,1]$ and orientation $(\theta,\phi)$.~\citep{Haeberlen:1976, Conturo:1996} Tensor-valued diffusion acquisition schemes have since been used to further investigate signal representations~\citep{Lasic:2014, Westin:2016, Cottaar:2020} and models.~\citep{Lampinen_CODIVIDE:2017, Coelho:2019, Coelho:2019_MICCAI, Reisert:2019, Lampinen:2020, Reymbaut_arxiv_Magic_DIAMOND:2020}

While inversion of the diffusion-relaxation NMR signal is already common practice in the porous media field,\citep{Prange:2009, Galvosas:2010, Bernin_Topgaard:2013, Song_book_chapter:2017} the advent of `tensor-valued' diffusion-relaxation correlation measurements has resulted in the development of non-parametric Monte-Carlo signal inversion algorithms of the diffusion~\citep{deAlmeidaMartins_Topgaard:2016} and diffusion-$T_1$-$T_2$~\citep{deAlmeidaMartins_Topgaard:2018} MRI signals in porous media, and of the diffusion~\citep{Topgaard:2019} and diffusion-$T_2$~\citep{deAlmeidaMartins:2020} MRI signals in the \textit{in vivo} brain.
Although noise-sensitive,~\citep{Reymbaut_accuracy_precision:2020} these algorithms do not rely on any compartmental/functional assumption regarding the voxel content, nor on constraints regarding data compression~\citep{Venkataramanan:2002} or the density of the acquisition sampling scheme.~\citep{Benjamini:2016, Kim:2017, Benjamini:2018, Benjamini:2020, Kim:2020} They also do not consider any regularization to guide the search for a suitable solution to the inverse problem.~\citep{Provencher:1982,Kroeker:1986,Whittall:1989,Mitchell:2012}
Enhanced by methods aiming to visualize and quantify fiber-specific properties, even in fiber-crossing areas of the white matter (60 to 90\% of voxels in a typical whole-brain imaging experiment~\citep{Jeurissen:2013}), Monte-Carlo signal inversions have been shown to yield critical sensitivity and specificity to fiber-specific $T_2$ values.~\citep{deAlmeidaMartins_ODF:2020, deAlmeidaMartins_thesis:2020, Reymbaut_arxiv_MC_DPC:2020} However, this work has yet to be extended to fiber-specific $\mathit{T}_1$-values, which are of particular interest to evaluate changes in bundle-specific myelin contents,~\citep{Liu:2019} relevant to the study of neurodevelopment, plasticity, aging and neurological disorders.~\citep{van_den_Heuvel:2010,Caeyenberghs:2016,Mancini:2018} Indeed, $T_1$ contrast is sensitive to myelin,~\citep{Mottershead:2003, Bjarnason:2005, deSantis:2014, Lutti:2014, Stuber:2014} like many other contrasts~\citep{Laule:2007, Campbell:2018} such as $T_2$,~\citep{Mackay:1994, Beaulieu:1998, Gareau:2000, Webb:2003, Stanisz:2004, Laule:2006, Laule:2008} $T_2^*$~\citep{Hwang:2010,Lee:2012,Sati:2013} and magnetization transfer.~\citep{Gareau:2000, Schmierer:2008} Importantly, note that none of these contrasts are `specific' to myelin.~\citep{Campbell:2018}

In this work, we leverage non-parametric distributions $\mathcal{P}(\mathbf{D},R_1)$ of diffusion tensors $\mathbf{D}$ and longitudinal relaxation rates $R_1 = 1/T_1$ obtained \textit{via} Monte-Carlo inversion of a diffusion-$T_1$ weighted \textit{in vivo} human-brain dataset to resolve sub-voxel diffusion-$R_1$ components. We first estimate parameter maps of these distributions' statistical descriptors and extract orientation-resolved $T_1$ values within the pool of highly anisotropic components output by the Monte-Carlo inversion algorithm. These $T_1$ values are then color-mapped onto non-parametric orientation distribution functions (ODFs)~\citep{deAlmeidaMartins_ODF:2020, deAlmeidaMartins_thesis:2020} and quantified in terms of median value and precision using orientational clusters obtained \textit{via} Monte-Carlo density-peak clustering (MC-DPC).~\citep{Reymbaut_arxiv_MC_DPC:2020} In particular, these novel tools enable to identify significant differences with respect to $T_1$ relaxation between major white-matter bundles. After describing how our \textit{in vivo} human-brain data was acquired in Section~\ref{Sec_in_vivo_data}, we lay down the theory underlying the Monte-Carlo signal inversion algorithm, the statistical descriptors of $\mathcal{P}(\mathbf{D},R_1)$, and our ODFs in Section~\ref{Sec_Monte_Carlo}, and detail the MC-DPC procedure in Section~\ref{Sec_MC_DPC}. We then present our results in Section~\ref{Sec_Results} and discuss them in Section~\ref{Sec_Discussion}, before concluding in Section~\ref{Sec_Conclusions}.

We emphasize that this work is merely a proof of concept for non-parametric fiber-specific $T_1$ relaxometry. Firstly, vast improvements could be brought to our experimental setup, as detailed in Sections~\ref{Sec_in_vivo_data} and \ref{Sec_Discussion}. Second, better matching between the output of our Monte-Carlo framework and plausible white-matter tracts would be yielded upon integrating tractography~\citep{Mori:1999,Basser:2000,Morris:2008,Reisert:2011,Christiaens:2015,Neher:2017,Konopleva:2018,Poulin:2019} into our analysis pipeline.

\section{Methods}
\label{Sec_Methods}

\subsection{In vivo human-brain data}
\label{Sec_in_vivo_data}

Data collection was approved by the Spectrum Medical Imaging local ethics committee. A healthy volunteer was scanned on a 3T GE 750w equipped with a 32-channel receiver head and neck GEM coils (only 12-16 channels used for head) using a prototype GE multidimensional diffusion (MDD) spin-echo sequence with EPI readout, echo time $\tau_\mathrm{E} =120$ ms, FOV=240x240x12 mm$^3$, voxel-size=3x3x3 mm$^3$, fat-saturation (fat-sat) pulses,~\cite{Haase:1985} and ASSET acceleration factor=2, customized for tensor-valued diffusion encoding~\citep{Lasic:2014,Szczepankiewicz_DIVIDE:2019} and variable repetition time $\tau_\mathrm{R}$. Tensor-valued diffusion encoding was performed with numerically optimized~\citep{Sjolund:2015} Maxwell-compensated~\citep{Szczepankiewicz_Maxwell:2019} waveforms. We also attempted to match their frequency contents.~\citep{Lundell:2019} The same tensor-valued diffusion-weighted sequence was repeated for $\tau_\mathrm{R} =1$, 2 and 5 s. The dimensions of the resulting 20-minute 363-point acquisition scheme, shown in Figure~\ref{Figure_acq}, match those of $\mathcal{P}(\mathbf{D},R_1)$. The signal-to-noise ratio of this dataset was estimated across voxels of the corona radiata by computing the mean-to-standard-deviation ratio of the spherically encoded diffusion signal at $b=0.1$ ms/{\textmu}m$^2$ (see Supplemental Material of Ref.~\onlinecite{Szczepankiewicz_DIVIDE:2019}). It equals 20 at $\tau_\mathrm{R}=1$ s and 40 at $\tau_\mathrm{R}=5$ s.

As indicated by the aforementioned FOV and voxel size, only four axial slices were acquired so that to limit the acquisition time. While an inversion-recovery slice-shuffling sequence could drastically reduce acquisition time,~\citep{Hutter:2018, Park_ISMRM:2018} our prototype sequence is currently limited to sequential slices.

\begin{figure*}[ht!]
\begin{center}
\includegraphics[width=30pc]{1_Fig_acq.pdf}
\caption{Visualization of our acquisition scheme. (A) Spin-echo sequence with EPI readout customized for tensor-valued diffusion encoding and variable repetition time $\tau_\mathrm{R}$. (B) 5D grid-like acquisition scheme where black points indicate the acquisition points in the 3D subspace of repetition time $\tau_\mathrm{R}$, b-tensor size $b$ and b-tensor shape $b_\Delta \in [-0.5, 1]$. The number of b-tensor orientations $(\Theta,\Phi)$ for each point is illustrated by the projected contours. (C) Acquisition parameters as a function of sorted acquisition point index $n_\mathrm{acq}$.}
\label{Figure_acq}
\end{center}
\end{figure*}

\subsection{Non-parametric Monte-Carlo inversion}
\label{Sec_Monte_Carlo}

\subsubsection{Signal fitting and bootstrapping}
\label{Sec_MC_algorithm}

We used a modified version of the 5D Monte-Carlo inversion algorithm found in Ref.~\onlinecite{deAlmeidaMartins:2020} (pioneered by Ref.~\onlinecite{Prange:2009}) to analyze the diffusion-$R_1$ dataset described in Section~\ref{Sec_in_vivo_data}. Let us consider axisymmetric diffusion tensors, parametrized by their axial diffusivity $\mathit{D}_\parallel$, radial diffusivity $\mathit{D}_\perp$ and orientation $(\theta, \phi)$. An alternative parametrization includes the isotropic diffusivity $\mathit{D}_\mathrm{iso} = (\mathit{D}_\parallel + 2\mathit{D}_\perp)/3$ and normalized anisotropy $\mathit{D}_\Delta = (\mathit{D}_\parallel - \mathit{D}_\perp)/(\mathit{D}_\parallel + 2\mathit{D}_\perp) \in [-0.5, 1]$.~\citep{Haeberlen:1976, Conturo:1996, Eriksson:2015} Our Monte-Carlo inversion technique retrieves non-parametric 5D intra-voxel distributions $\mathcal{P}(\mathbf{D},\mathit{R}_1) \equiv \mathcal{P}(\mathit{D}_\parallel,\mathit{D}_\perp,\theta,\phi, \mathit{R}_1)$ by fitting a set of diffusion-$T_1$ weighted signals with a finite weighted sum of $\mathit{N}$ components $(\mathbf{D}_\mathit{n}, \mathit{R}_{1,\mathit{n}})\equiv (\mathit{D}_{\parallel,\mathit{n}}, \mathit{D}_{\perp,\mathit{n}} , \theta_\mathit{n}, \phi_\mathit{n}, \mathit{R}_{1,\mathit{n}})$, with $1\leq n\leq N=50$. Given that the $\mathit{T}_1$-weighting of the dataset detailed in Section~\ref{Sec_in_vivo_data} is provided through a spin-echo sequence with constant echo time $\tau_\mathrm{E}$ and variable repetition time $\tau_\mathrm{R}$, the inversion algorithm inverts the following discretized signal equation~\citep{Perman:1984}
\begin{align}
\mathcal{S}_\mathit{m} = & \sum_{\mathit{n}=1}^\mathit{N} \mathit{w}_\mathit{n}\,\exp(-\mathbf{b}_\mathit{m}:\mathbf{D}_\mathit{n}) \nonumber \\
 & \times \left[1 -2\exp([\tau_{\mathrm{E}}/2-\tau_{\mathrm{R},\mathit{m}}]\mathit{R}_{1,\mathit{n}}) + \exp(-\tau_{\mathrm{R},\mathit{m}}\mathit{R}_{1,\mathit{n}})\right]\, ,
\label{Eq_signal_discretized}
\end{align}
where $\mathcal{S}_\mathit{m}$ is the $m^\text{th}$ acquired signal, $\mathit{w}_\mathit{n}$ is the weight of the $n^\text{th}$ component, $\mathbf{b}$ is the diffusion-encoding tensor (b-tensor)~\citep{Mattiello:1994, Mattiello:1997} from tensor-valued diffusion encoding,~\citep{Eriksson:2013,Westin:2014,Eriksson:2015,Westin:2016,Topgaard:2017,Topgaard_dim_rand_walks:2019} and ":" is the Frobenius inner product. The weights $w_n$ are normalized so that $\sum_{\mathit{n}=1}^\mathit{N} \mathit{w}_\mathit{n} = \mathcal{S}_0 = \mathcal{S}(\mathbf{b}=\mathbf{0},\mathrm{\tau}_\mathrm{R}\to +\infty)$. For axisymmetric b-tensors, the Frobenius inner product writes~\citep{Eriksson:2015} $\mathbf{b}:\mathbf{D} = bD_\mathrm{iso}[1 + 2b_\Delta D_\Delta P_2(\cos\beta)]$, where $\mathit{P}_2(\mathit{x}) = (3\mathit{x}^2-1)/2$ is the second Legendre polynomial and $\cos\beta = \cos\Theta \cos\theta + \sin\Theta\sin\theta\cos(\Phi - \phi)$ is the cosine of the shortest angle $\beta$ between the main axis $(\Theta,\Phi)$ of $\mathbf{b}$ and the main axis $(\theta,\phi)$ of $\mathbf{D}$.

A short-hand notation of Equation~\ref{Eq_signal_discretized} reads
\begin{equation}
\mathbf{S} = \mathbf{K}\mathbf{w}\, ,
\end{equation}
where $\mathbf{S}$ is the column vector containing the acquired signals $\mathcal{S}_m$, $\mathbf{K}$ is the inversion kernel matrix containing the signal decays and $\mathbf{w}$ is the column vector containing the weights $\mathit{w}_n$ of the components $(\mathbf{D}_n, \mathit{R}_{1,n})$.
The Monte-Carlo inversion algorithm randomly samples such components within the following ranges, $D_\parallel,D_\perp \in [5\times 10^{-3}, 5]\;\text{\textmu}\mathrm{m}^2/\mathrm{ms}$, $\cos\theta\in[0,1[$, $\phi\in[0,2\pi[$ and $R_1\in[0.1,2]\;\mathrm{s}^{-1}$, and estimates the associated vector $\mathbf{w}$ quantifying the components' propensity to fit the acquired signals \textit{via} non-negative least-squares fitting:~\citep{Lawson_book:1974,Whittall:1989,English:1991,Venkataramanan:2002,Mitchell:2012}
\begin{equation}
\mathbf{w} = \underset{\mathbf{w}^\prime\geq 0}{\mathrm{argmin}}\;\Vert\mathbf{S}-\mathbf{K}\mathbf{w}^\prime\Vert_2^2\, ,
\end{equation}
where $\Vert\cdot\Vert_2$ denotes the L$_2$ norm. This process is repeated iteratively following a quasi-genetic filtering detailed in Refs.~\onlinecite{deAlmeidaMartins_Topgaard:2016, deAlmeidaMartins_Topgaard:2018, Topgaard:2019, deAlmeidaMartins:2020}.
Embracing the inherent ill-conditioning of Laplace inversion problems, we performed bootstrapping with replacement~\citep{de_Kort:2014} on the data and estimated for each voxel an ensemble of $\mathit{N}_\mathrm{b}=96$ plausible sets of components, also called "bootstrap solutions", each denoted by $\{(\mathit{D}_{\parallel,\mathit{n}}, \mathit{D}_{\perp,\mathit{n}} , \theta_\mathit{n}, \phi_\mathit{n}, \mathit{R}_{1,\mathit{n}}, \mathit{w}_\mathit{n})\}_{1\leq \mathit{n}\leq \mathit{N}=20}$. We then computed statistical descriptors of $\mathcal{P}(\mathbf{D},\mathit{R}_1)$ for each bootstrap solution and calculated the median of each statistical descriptor across bootstrap solutions (see Section~\ref{Sec_statistical_descriptors_binning}).

\subsubsection{Statistical descriptors and binning}
\label{Sec_statistical_descriptors_binning}

The final solution of the Monte-Carlo inversion algorithm, $\mathcal{P}(\mathbf{D},\mathit{R}_1)$, can be understood as the median of the solutions obtained for each bootstrap solution, $\mathcal{P}_{\mathit{n}_\mathrm{b}}(\mathbf{D},\mathit{R}_1)$, with $1\leq \mathit{n}_\mathrm{b} \leq \mathit{N}_\mathrm{b}=96$. Following previous works,~\citep{Reymbaut_accuracy_precision:2020,deAlmeidaMartins:2020} we quantified the main features of this final solution by computing the median across bootstrap solutions of means $\mathrm{Med}_{(n_\mathrm{b})}\,(\mathrm{E}[\chi]_{n_\mathrm{b}})$, variances $\mathrm{Med}_{(n_\mathrm{b})}\,(\mathrm{V}[\chi]_{n_\mathrm{b}})$ and covariances $\mathrm{Med}_{(n_\mathrm{b})}\,(\mathrm{C}[\chi,\chi^\prime]_{n_\mathrm{b}})$ of the per-bootstrap isotropic diffusivity, squared normalized anisotropy, and longitudinal relaxation rate $\chi,\chi^\prime=\mathit{D}_\mathrm{iso}, D_\Delta^2, R_1$, respectively. Here, the median operator $\mathrm{Med}_{(n_\mathrm{b})}(\,\cdot\,)$ acts across bootstrap solutions and $\mathrm{E}[\,\cdot\,]_{\mathit{n}_\mathrm{b}}$, $\mathrm{V}[\,\cdot\,]_{\mathit{n}_\mathrm{b}}$ and $\mathrm{C}[\,\cdot,\cdot\,]_{\mathit{n}_\mathrm{b}}$ denote the per-bootstrap average, variance and covariance over the diffusion-relaxation components forming the bootstrap solution $\mathit{n}_\mathrm{b}$, respectively. For simplicity, we omit the explicit mention of the median operator when addressing a statistical descriptor, thereby writing averages, variances and covariances as $\mathrm{E}[\chi]$, $\mathrm{V}[\chi]$ and $\mathrm{C}[\chi,\chi^\prime]$ respectively.

By design, the Monte-Carlo inversion algorithm progressively builds up the sought-for intra-voxel distribution $\mathcal{P}(\mathbf{D},\mathit{R}_1)$ as a non-parametric discrete weighted sum of components. This implies that tissue-specific statistical descriptors can be extracted by subdividing the 5D configuration space of $\mathcal{P}(\mathbf{D},\mathit{R}_1)$ into multiple bins. For instance, the "thin", "thick" and "big" bins introduced in Refs.~\onlinecite{Topgaard:2019, deAlmeidaMartins:2020} aim to isolate the signal contributions from white matter, grey matter and cerebrospinal fluid, respectively. In this work, the boundaries of these bins, illustrated in the panels C, D and E of Figure~\ref{Figure_fit_distributions}, were defined as follows:
\begin{itemize}
\item "big" bin within $D_\mathrm{iso} \in [2, 10] \;\text{\textmu}\mathrm{m}^2/\mathrm{ms}$, $D_\parallel/D_\perp \in [0.01, 1000]$ and $R_1 \in [0.01, 10] \;\mathrm{s}^{-1}$.
\item "thick" bin within $D_\mathrm{iso} \in [0.1, 2] \;\text{\textmu}\mathrm{m}^2/\mathrm{ms}$, $D_\parallel/D_\perp \in [0.01, 4]$ and $R_1 \in [0.01, 10] \;\mathrm{s}^{-1}$.
\item "thin" bin within $D_\mathrm{iso} \in [0.1, 2] \;\text{\textmu}\mathrm{m}^2/\mathrm{ms}$, $D_\parallel/D_\perp \in [4, 1000]$ and $R_1 \in [0.01, 10] \;\mathrm{s}^{-1}$.
\end{itemize}
As such, the "big" bin captures highly diffusive components, the "thick" bin captures components that are not highly diffusive nor highly anisotropic, and the "thin" bin captures components that are highly anisotropic. Note that the above bin boundaries extend far beyond the inversion boundaries listed in Section~\ref{Sec_MC_algorithm}, so that to produce aesthetically pleasing figures such as the panels C, D and E of Figure~\ref{Figure_fit_distributions}. Bin-specific statistical descriptors were estimated following the aforementioned process for the retrieved components that specifically fall into each bin.

\subsubsection{Orientation distribution functions (ODFs)}
\label{Sec_ODFs}

Orientation distribution functions (ODFs) can be defined from the thin-bin components output by the Monte-Carlo signal inversion of Section~\ref{Sec_MC_algorithm} using the procedure detailed in Refs.~\onlinecite{deAlmeidaMartins_ODF:2020, deAlmeidaMartins_thesis:2020}. Briefly, per-bootstrap ODFs were generated by convolving the discrete ensemble of weights $w_i$ and unit-orientations $(\theta_i, \phi_i)$ of the components belonging to the thin bin in each bootstrap solution $n_\mathrm{b}$, $\{(w_\mathit{i}, \theta_i, \phi_i)\}_{\mathit{n}_\mathrm{b},\;\mathit{i}\in \{\text{thin bin}\}}$, with a Watson kernel. The purpose of this kernel is to smoothly map the discrete set of components onto the nearest nodes of a dense spherical mesh $\{(\theta_{\text{mesh}},\phi_{\text{mesh}})\}$. A final voxel-wise ODF $\mathit{P}(\theta_{\text{mesh}},\phi_{\text{mesh}})$ was calculated as the median of the per-bootstrap ODFs. In this work, we considered a $1000$-point uniform spherical mesh and set the concentration parameter of the Watson kernel to $\kappa=14.9$, following the rationale detailed in Ref.~\onlinecite{deAlmeidaMartins_ODF:2020}.

Following a similar procedure, the discrete set of diffusion-relaxation metrics $\{(\mathit{D}_{\parallel,\mathit{i}}, \mathit{D}_{\perp,\mathit{i}}, \mathit{R}_{1,\mathit{i}})\}_{\mathit{n}_\mathrm{b},\;\mathit{i}\in \{\text{thin bin}\}}$ can also be mapped onto this spherical mesh,~\citep{deAlmeidaMartins_ODF:2020, deAlmeidaMartins_thesis:2020} yielding the orientation-specific diffusion-relaxation measures $\hat{\mathrm{E}}[\chi]_{\mathit{n}_\mathrm{b}}(\theta_\text{mesh},\phi_\text{mesh})$, with $\chi \equiv \mathit{D}_\mathrm{iso}, \mathit{D}_\Delta^2, \mathit{R}_1$. In addition to coloring ODFs according to the $(\theta_{\text{mesh}},\phi_{\text{mesh}})$ local orientation, this mapping allows to color ODFs according to the local value $\mathrm{Med}_{(n_\mathrm{b})}\left(\hat{\mathrm{E}}[\chi]_{\mathit{n}_\mathrm{b}}(\theta_\text{mesh},\phi_\text{mesh})\right)$, thereby improving the visualization of orientation-specific diffusion-relaxation quantities. For simplicity, the short-hand notation "$\hat{\mathrm{E}}[\chi]$" is now retained instead for $\mathrm{Med}_{(n_\mathrm{b})}\left(\hat{\mathrm{E}}[\chi]_{\mathit{n}_\mathrm{b}}(\theta_\text{mesh},\phi_\text{mesh})\right)$.

\subsection{Monte-Carlo density-peak clustering (MC-DPC)}
\label{Sec_MC_DPC}

The Monte-Carlo signal inversion algorithm of Section~\ref{Sec_MC_algorithm} can be combined with density-peak clustering (DPC)~\citep{Rodriguez_Laio:2014} according to the work presented in Ref.~\onlinecite{Reymbaut_arxiv_MC_DPC:2020}. This combination, called "Monte-Carlo density-peak clustering" (MC-DPC), enables to quantify the median value and precision of orientation-resolved means of $\chi=\mathit{D}_\mathrm{iso}, D_\Delta^2, R_1, T_1=1/R_1$ across bootstrap solutions. We used MC-DPC to detect statistically significant differences between sub-voxel fiber populations robustly assigned to major fiber bundles \textit{a posteriori}.

Firstly, MC-DPC gathers the ensemble $\{ \mathcal{E}^\mathrm{thin}_{n_\mathrm{b}} \}_{1\leq n_\mathrm{b}\leq N_\mathrm{b}}$ of all per-bootstrap thin-bin solution sets $\mathcal{E}^\mathrm{thin}_{n_\mathrm{b}} = \{(\mathit{D}_{\parallel,\mathit{i}}, \mathit{D}_{\perp,\mathit{i}} , \theta_\mathit{i}, \phi_\mathit{i}, \mathit{R}_{1,\mathit{i}}, \mathit{w}_\mathit{i})\}_{n_\mathrm{b},\;\mathit{i}\in \{\text{thin bin}\}}$ and delineates $N_\mathrm{c}$ clusters in its orientation subspace using DPC with data-point density and outlier detection altered to account for the weights $w_i$ of the retrieved thin-bin components.~\citep{Reymbaut_arxiv_MC_DPC:2020} An initial number of clusters $N_\mathrm{c}$ is automatically set by the number of voxel-wise ODF peaks, but may be reduced by MC-DPC following a filtering approach detailed in Ref.~\onlinecite{Reymbaut_arxiv_MC_DPC:2020}. Assuming that the estimated clusters, resulting from orientational aggregates of the all-bootstrap thin-bin solutions, can be interpreted as orientational regions of interest associated with sub-voxel fiber populations, MC-DPC then computes orientation-resolved statistics across bootstrap solutions. To do so, it separately classifies each per-bootstrap ensemble of thin-bin solutions $\mathcal{E}^\mathrm{thin}_{n_\mathrm{b}}$ into $N_\mathrm{c}$ ensembles $\mathcal{E}^\mathrm{thin}_{n_\mathrm{b},n_\mathrm{c}}$ (with $1\leq n_\mathrm{c}\leq N_\mathrm{c}$), each containing the thin-bin solutions of bootstrap solution $n_\mathrm{b}$ that belong to an estimated cluster $n_\mathrm{c}$. It then averages the properties of the solutions within each ensemble $\mathcal{E}^\mathrm{thin}_{n_\mathrm{b},n_\mathrm{c}}$ independently, yielding the orientation-resolved means
\begin{equation}
\mathring{\mathrm{E}}[\chi]_{n_\mathrm{b},n_\mathrm{c}} = \frac{\sum_{k\in  \mathcal{E}^\mathrm{thin}_{n_\mathrm{b},n_\mathrm{c}}} w_k\,\chi_k}{\sum_{k\in \mathcal{E}^\mathrm{thin}_{n_\mathrm{b},n_\mathrm{c}}} w_k} \, ,
\label{Eq_orientation_resolved_mean}
\end{equation}
with $\chi\equiv x,y,z, \mathit{D}_\mathrm{iso}, \mathit{D}_\Delta^2, \mathit{R}_1, \mathit{T}_1=1/\mathit{R}_1$, where $(x,y,z)$ are the Cartesian coordinates of a component's orientation $(\theta,\phi)$. The short-hand notation "$\mathring{\mathrm{E}}[\chi]$" will be used for simplicity to describe the collection of orientation-resolved means $\mathring{\mathrm{E}}[\chi]_{n_\mathrm{b},n_\mathrm{c}}$ originating from all bootstrap solutions $n_\mathrm{b}$ and all clusters $n_\mathrm{c}$. The orientation-resolved means of $T_1$ and $R_1$ in Equation~\ref{Eq_orientation_resolved_mean} are computed separately, as both quantities are commonly found in the MRI literature and $\mathring{\mathrm{E}}[T_1]$ does not equal $1/\mathring{\mathrm{E}}[R_1]$. Finally, we extracted the median and interquartile range of our various orientation-resolved means $\mathring{\mathrm{E}}[\chi]$ across bootstrap solutions.

Non-parametric Mann-Whitney $U$-tests were used to assess whether or not two orientation-resolved means $\mathring{\mathrm{E}}[\chi]$ characterizing distinct white-matter fiber bundles are sampled from identically shaped non-median-shifted continuous distributions (null hypothesis $\mathcal{H}_0$). The $p$-values resulting from these tests inform on the acceptance or rejection of $\mathcal{H}_0$ at a certain significance level. In particular, rejections of $\mathcal{H}_0$ at $p < 0.01$, $0.01\leq p < 0.05$ and $0.05\leq p < 0.1$ were used as proxies for detecting significant differences between two bundles with respect to a given type of orientation-resolved mean.

\section{Results}
\label{Sec_Results}

Figure~\ref{Figure_fit_distributions} presents the fitted signals and the intra-voxel distributions $\mathcal{P}(\mathbf{D},R_1)$ estimated by the Monte-Carlo signal inversion algorithm of Section~\ref{Sec_MC_algorithm} in typical voxels associated with cerebrospinal fluid (CSF) in the ventricles, cortical grey matter (GM) and white matter (WM) in the corpus callosum. Figure~\ref{Figure_maps} displays typical axial maps of the statistical descriptors and bin-specific statistical descriptors described in Section~\ref{Sec_statistical_descriptors_binning}. Figure~\ref{Figure_ODFs} shows orientation-colored and $\hat{\mathrm{E}}[R_1]$-colored ODFs (see Section~\ref{Sec_ODFs}) in a typical axial slice. Figures~\ref{Figure_CC_CING} and \ref{Figure_CS} investigate possible microstructural differences between sub-voxel fiber populations leveraging MC-DPC (see Section~\ref{Sec_MC_DPC}) in regions of interest that target specific fiber crossings. While Figure~\ref{Figure_CC_CING} focuses on the crossing between the corpus callosum (CC) and the cingulum (CING), Figure~\ref{Figure_CS} focuses on the crossing between the corpus callosum (CC), the arcuate fasciculus (AF) and the corticospinal tract (CST) in the posterior corona radiata.

\begin{figure*}[ht!]
\begin{center}
\includegraphics[width=30pc]{2_Fig_fit_distributions.pdf}
\caption{Monte-Carlo fitted signal and retrieved 5D distributions $\mathcal{P}(\mathbf{D},R_1)$ in typical voxels. (A) $\mathcal{S}_0$ map estimated by the Monte-Carlo inversion. The colored squares delineate typical CSF (blue), GM (green) and WM (red) voxels. (B) Normalized signal $\tilde{S} = \mathcal{S}/\mathrm{max}(\mathcal{S})$ measured (black points) and fitted (colored points) in the archetypal voxels of panel A as a function of the sorted acquisition point index $n_\mathrm{acq}$ of Figure~\ref{Figure_acq}. (C,D,E) Non-parametric distributions $\mathcal{P}(\mathbf{D},R_1)$ estimated within the archetypal voxels of panel A and reported as scatter plots in a 3D space of the logarithms of the longitudinal relaxation rate $R_1$, isotropic diffusivity $D_\mathrm{iso}$, and axial-radial diffusivity ratio $D_\parallel/D_\perp$. Diffusion-orientations $(\theta, \phi)$ are color-coded according to $[\mathrm{red,green,blue}] = [\sin\theta\cos\phi, \sin\theta\sin\phi, \cos\theta]\times \vert D_\parallel - D_\perp\vert/\mathrm{max}(D_\parallel,D_\perp)$. Symbol area is proportional to the statistical weight $w_n/\mathcal{S}_0$ of the corresponding component $n$. The contour lines on the sides of the plots represent projections of the 5D distributions $\mathcal{P}(\mathbf{D},R_1)$ onto the respective 2D planes. The "thin", "thick" and "big" bins defined in Section~\ref{Sec_statistical_descriptors_binning} are illustrated in the panels where they are most relevant.}
\label{Figure_fit_distributions}
\end{center}
\end{figure*}

\begin{figure*}[ht!]
\begin{center}
\includegraphics[width=30pc]{3_Fig_maps.pdf}
\caption{Typical axial maps of the statistical descriptors described in Section~\ref{Sec_statistical_descriptors_binning}. The various colormaps are bounded by the following values: $\mathcal{S}_0\in[0, \mathrm{max}(\mathcal{S}_0)]$, $\mathrm{E}[D_\mathrm{iso}]\in[0, 3.5]\;\text{\textmu}\mathrm{m}^2/\mathrm{ms}$, $\mathrm{E}[D_\Delta^2]\in[0, 1]$, $\mathrm{E}[R_1]\in [0,0.8]\;\mathrm{s}^{-1}$, $\mathrm{V}[D_\mathrm{iso}]\in[0, 0.9]\;\text{\textmu}\mathrm{m}^4/\mathrm{ms}^2$, $\mathrm{V}[D_\Delta^2]\in[0, 0.1]$, $\mathrm{V}[R_1]\in[0, 0.13]\;\mathrm{s}^{-2}$, $\mathrm{C}[D_\mathrm{iso},D_\Delta^2]\in[-0.3, 0.3]\;\text{\textmu}\mathrm{m}^2/\mathrm{ms}$, $\mathrm{C}[D_\mathrm{iso},R_1]\in[-2.4, 2.4]\times 10^{-4}\;\text{\textmu}\mathrm{m}^2/\mathrm{ms}^2$ and $\mathrm{C}[D_\Delta^2,R_1]\in[-0.08, 0.08]\;\mathrm{s}^{-1}$. The bin-specific average intra-voxel orientation $\mathrm{E}[(\theta,\phi)]$ is color-coded for orientation according to $[\mathrm{red,green,blue}] = [\mathrm{E}[D_{xx}], \mathrm{E}[D_{yy}], \mathrm{E}[D_{zz}]]/\mathrm{max}(\mathrm{E}[D_{xx}], \mathrm{E}[D_{yy}], \mathrm{E}[D_{zz}])$, where the average diffusivities $\mathrm{E}[D_{ii}]$ are associated with the directions $i=x, y, z$ corresponding to the "left-right", "anterior-posterior" and "superior-inferior" directions, respectively. For a given bin, the intensity of the bin-specific maps equals the voxel-wise average fraction $f_\mathrm{bin}$ of components belonging to this bin. The segmentation map is colored according to $[\mathrm{red,green,blue}] = [f_\mathrm{thin}, f_\mathrm{thick}, f_\mathrm{big}]/\mathrm{max}(f_\mathrm{thin}, f_\mathrm{thick}, f_\mathrm{big})$.}
\label{Figure_maps}
\end{center}
\end{figure*}

\begin{figure*}[ht!]
\begin{center}
\includegraphics[width=30pc]{4_Fig_ODFs.pdf}
\caption{Axial greyscale maps of the fraction of non thin-bin components $1-f_\mathrm{thin}$ with superimposed ODFs colored by local orientation (with $x$, $y$ and $z$ corresponding to the "left-right", "anterior-posterior" and "superior-inferior" directions, respectively) and by $\hat{\mathrm{E}}[R_1]$ (see Section~\ref{Sec_ODFs}). The middle insets zoom on a voxel containing a fiber crossing between the corpus callosum (CC) and the cingulum (CING), and presents the estimated orientation-colored and $\hat{\mathrm{E}}[R_1]$-colored ODFs for this voxel. While differences in $\hat{\mathrm{E}}[R_1]$ seem to exist between CC and CING, such differences may also exist in the regions where the CST's pyramidal tracts are located (blue circles), as indicated by the greener $\hat{\mathrm{E}}[R_1]$-colored ODFs therein.}
\label{Figure_ODFs}
\end{center}
\end{figure*}

\begin{figure*}[ht!]
\begin{center}
\includegraphics[width=33pc]{5_Fig_CC_CING.pdf}
\caption{Boxplots of the medians of the orientation-resolved means $\mathring{\mathrm{E}}[\chi]$ obtained from MC-DPC (see Equation~\ref{Eq_orientation_resolved_mean}) within the regions of interest drawn in the left panel (white-lined boxes), showing two axial slices of the orientation-colored average fraction of thin-bin components $f_\mathrm{thin}$. For a given boxplot, the horizontal line and whiskers indicate the median and the range between the first and third quartiles of the medians of the orientation-resolved means $\mathring{\mathrm{E}}[\chi]$, respectively. The chosen regions of interest focus on crossing areas between the corpus callosum (CC) and the cingulum (CING). Each MC-DPC cluster (and associated orientation-resolved means) is robustly assigned to one of these bundles depending on whether its median orientation is closer to the $x$ "left-right" direction (CC) or to the $y$ "anterior-posterior" direction (CING). The asterisks report the results of non-parametric Mann-Whitney $U$-tests assessing whether or not two orientation-resolved means $\mathring{\mathrm{E}}[\chi]$ assigned to distinct bundles are sampled from identically shaped non-median-shifted continuous distributions (null hypothesis $\mathcal{H}_0$). The $p$-values resulting from these tests inform on the acceptance or rejection of $\mathcal{H}_0$ at a certain significance level: $0.05\leq p < 0.1$ (*) and $0.01\leq p < 0.05$ (**).}
\label{Figure_CC_CING}
\end{center}
\end{figure*}

\begin{figure*}[ht!]
\begin{center}
\includegraphics[width=33pc]{6_Fig_CS.pdf}
\caption{Equivalent of Figure~\ref{Figure_CC_CING} for regions of interest focusing on crossing areas between the corpus callosum (CC), the arcuate fasciculus (AF) and the corticospinal tract (CST) in the posterior corona radiata. Each MC-DPC cluster (and associated orientation-resolved means) is robustly assigned to one of these bundles depending on whether its median orientation is closer to the $x$ "left-right" direction (CC), to the $y$ "anterior-posterior" direction (AF), or to the $z$ "superior-inferior" direction (CST). The asterisks report the results of non-parametric Mann-Whitney $U$-tests assessing whether or not two orientation-resolved means $\mathring{\mathrm{E}}[\chi]$ assigned to distinct bundles are sampled from identically shaped non-median-shifted continuous distributions (null hypothesis $\mathcal{H}_0$). The $p$-values resulting from these tests inform on the acceptance or rejection of $\mathcal{H}_0$ at a certain significance level: $0.05\leq p < 0.1$ (*), $0.01\leq p < 0.05$ (**) and $p < 0.01$ (***).}
\label{Figure_CS}
\end{center}
\end{figure*}

\section{Discussion}
\label{Sec_Discussion}

Figure~\ref{Figure_fit_distributions} demonstrates that the Monte-Carlo signal inversion algorithm described in Section~\ref{Sec_MC_algorithm} adequately fits the measured signal in various voxels pertaining to WM, cortical GM and CSF in the ventricles. Additionally, the bins defined in Section~\ref{Sec_statistical_descriptors_binning} appear to capture these environments in accordance with their original design, \textit{i.e.} the thin, thick and big bins capture components typically associated with WM, GM and CSF, respectively.

Figure~\ref{Figure_maps} shows that the Monte-Carlo signal inversion algorithm can estimate maps of $\mathcal{P}(\mathbf{D},R_1)$'s statistical descriptors. In particular, it retrieves maps that are consistent with those thoroughly discussed in Ref.~\onlinecite{deAlmeidaMartins:2020}, namely the $\mathcal{S}_0$, $\mathrm{E}[D_\mathrm{iso}]$, $\mathrm{E}[D_\Delta^2]$, $\mathrm{V}[D_\mathrm{iso}]$, $\mathrm{V}[D_\Delta^2]$ and $\mathrm{C}[D_\mathrm{iso},D_\Delta^2]$ maps, along with the bin-specific $\mathrm{E}[D_\mathrm{iso}]$, $\mathrm{E}[D_\Delta^2]$ and average local orientation $\mathrm{E}[(\theta,\phi)]$ maps, and the bin-segmentation map. Let us thus mainly discuss the $R_1$-related maps. The $\mathrm{E}[R_1]$ map resembles an expected low-resolution $T_1$ map, \textit{i.e.} bright in WM, slightly darker in GM and very dark in CSF. The $\mathrm{V}[R_1]$ map resembles a noisier version of the $\mathrm{V}[D_\mathrm{iso}]$ map, as both give high values in mixed CSF-WM/GM voxels. The noise in $\mathrm{V}[R_1]$ could be reduced upon adding more repetition times in the acquisition scheme. The $\mathrm{C}[D_\mathrm{iso},R_1]$ map is negative at the interface between CSF and either WM or cortical GM. Indeed, upon entering CSF from WM/GM, $D_\mathrm{iso}$ increases rapidly and $T_1$ relaxation slows down (rapidly decreasing $R_1$). Finally, the $\mathrm{C}[D_\Delta^2,R_1]$ map exhibits no specific pattern. As for the bin-specific $\mathrm{E}[R_1]$ maps, they present a clear contrast between our tissue-specific bins, due to a fast $T_1$ relaxation in WM, a slower relaxation in GM and an even slower relaxation in CSF.

Figure~\ref{Figure_ODFs} features non-parametric ODFs capturing local orientations that are consistent with the known anatomy. Regarding $\hat{\mathrm{E}}[R_1]$-colored ODFs (see Section~\ref{Sec_ODFs}), they appear to change colors when approaching tissue interfaces with CSF. This gradual change in $\hat{\mathrm{E}}[R_1]$ may originate from exchange between tissues and CSF in these regions. Importantly, Figure~\ref{Figure_ODFs} shows that potential differences in $T_1$ relaxation may exist between major fiber bundles, namely CC and CING (insets in Figure~\ref{Figure_ODFs}), and CC, AF and CST (as indicated by greener $\hat{\mathrm{E}}[R_1]$-colored ODFs where the CST's pyramidal tracts are located).

These potential microstructural differences are quantified in Figures~\ref{Figure_CC_CING} and \ref{Figure_CS}. Let us focus on relaxation-based differences. Figure~\ref{Figure_CC_CING} shows that CC and CING exhibit significant differences in $\mathring{\mathrm{E}}[R_1]$ that are qualitatively consistent with those found in Refs.~\onlinecite{deSantis_T1:2016, Andrews_ISMRM:2019}, \textit{i.e.} $R_1$ tends to be lower in CING compared to CC. As for Figure~\ref{Figure_CS}, it shows that CST features significant differences in $\mathring{\mathrm{E}}[R_1]$ and $\mathring{\mathrm{E}}[T_1]$ with CC and AF (but no statistically significant differences between CC and AF). These differences are qualitatively consistent with those identified for CST in Ref.~\onlinecite{deSantis_T1:2016}, \textit{i.e.} $T_1$ tends to be higher in CST compared to CC and AF.

Quantitatively, the $T_1$ values estimated by $\mathring{\mathrm{E}}[T_1]$ in Figures~\ref{Figure_CC_CING} and \ref{Figure_CS} (around 1.5 to 2 seconds) are overestimated compared to those of Ref.~\onlinecite{deSantis_T1:2016} (around 0.9 to 1 second) and Ref.~\onlinecite{Andrews_ISMRM:2019} (around 0.7 second). This discrepancy can be explained by either or both of the following factors. Firstly, the acquisition scheme described in Section~\ref{Sec_in_vivo_data} does not does not maximize the amount of diffusion-relaxation correlations built into the inversion kernel of Equation~\ref{Eq_signal_discretized}, because the same diffusion-weighting block was repeated for each acquired repetition time. Similar problems have been suggested to lead to a loss of accuracy for the Monte-Carlo signal inversion.~\citep{Reymbaut_accuracy_precision:2020} Second, the use of saturation recovery with a spoiled spin echo for $T_1$ encoding is very sensitive to flip-angle inaccuracies caused by both $B_1^+$ inhomogeneity across the subject and slice profile imperfections. Saturation-recovery based $T_1$ mapping is also sensitive to magnetization-transfer effects, especially in the present setup comprising an additional refocusing pulse as well as a fat-saturation pulse.~\citep{Wolff_Balaban:1989, Teixeira:2019, Teixeira:2020} These technical limitations should be mitigated upon developing a sequence that includes inversion preparation for enhanced $T_1$ sensitivity and slice shuffling for optimized time efficiency.~\citep{Hutter:2018, Park_ISMRM:2018}

\newpage

\section{Conclusions}
\label{Sec_Conclusions}

Diffusion-$T_1$ weighted datasets incorporating multiple b-tensor shapes can be inverted to obtain non-parametric distributions $\mathcal{P}(\mathbf{D},R_1)$ of diffusion tensors and longitudinal relaxation rates using the Monte-Carlo signal inversion algorithm. The main features of the retrieved distributions can be visualized as maps and bin-specific maps of statistical descriptors related to means, variances and covariances of diffusion-relaxation properties. In particular, the bin-specific $\mathrm{E}[R_1]$ maps exhibit the expected $R_1$ contrast between white matter, grey matter and CSF. Further insight into white-matter microstructure is provided by the "thin bin", which isolates highly anisotropic components that should report on white-matter tissues. From these thin-bin components, visualization of fiber-specific information is improved upon defining orientation distribution functions (ODFs) that can be color-mapped with respect to local orientation or diffusion-relaxation features.~\citep{deAlmeidaMartins_ODF:2020} While $\hat{\mathrm{E}}[R_1]$-colored ODFs hint at possible differences between fiber bundles, Monte-Carlo density-peak clustering (MC-DPC) enables their quantification in terms of fiber-specific diffusion-relaxation measures.~\citep{Reymbaut_arxiv_MC_DPC:2020}

Importantly, significant differences with respect to longitudinal relaxation are detected between the corpus callosum and the cingulum, and between the corticospinal tract and the corpus callosum and arcuate fasciculus. These differences, qualitatively consistent with those found in previous works,~\citep{deSantis_T1:2016, Andrews_ISMRM:2019} offer a proof of concept for the potential of our Monte-Carlo framework in terms of non-parametric fiber-specific $T_1$ relaxometry. Such approach would be practical in identifying differences in $\mathit{T}_1$ between distinct sub-voxel fiber populations, characterizing developmental or pathological changes in $\mathit{T}_1$ within a given sub-voxel fiber population, and measuring the angular dependence of longitudinal relaxation times in white matter with respect to the main MRI magnetic field $\mathbf{B}_0$.~\citep{Henkelman:1994,Knight:2018} Moreover, non-parametric fiber-specific $T_1$ relaxometry would be highly relevant to microstructure-informed tractography.~\citep{Daducci:2016, Girard:2017, Barakovic_thesis:2019}

Nevertheless, our work can still be improved in two main ways. Firstly, the acquisition scheme could be optimized in terms of speed~\citep{Hutter:2018, Park_ISMRM:2018} and sensitivity.~\citep{Song:2005, Bates_ISMRM:2019, Song:2020} Second, MC-DPC could be combined with tractography to better assign MC-DPC's output clusters to their corresponding fiber bundles. These ideas will be explored in future investigations.

\newpage
\clearpage
\newpage

\section*{Acknowledgments}
This work was financially supported by the Swedish Foundation for Strategic Research (ITM17-0267) and the Swedish Research Council (2018-03697). D. Topgaard owns shares in Random Walk Imaging AB (Lund, Sweden, \href{http://www.rwi.se/}{http://www.rwi.se/}), holding patents related to the described methods.