Skip to content

HierarchicalProblem

Joint Bayesian inference over shared hyperparameters and group-level parameters across J nominally identical structures (or any grouped data).

hproblem = mc.HierarchicalProblem(
    hyperprior=hyperprior,          # log p(phi)
    group_prior=group_prior,        # log p(theta_j | phi)
    group_likelihoods=likelihoods,  # list of likelihoods, one per group
    n_hyper=2,
    n_group=1,
    param_names_hyper=["mu_k", "log_sigma_k"],
    param_names_group=["k"],
)

# Compatible with all samplers — run as a normal Problem
result = mc.DRAM(n_samples=20_000, initial_cov=np.eye(hproblem.n_params)).run(
    hproblem,
    x0=hproblem.default_x0(phi0=[10.0, 0.0], group_x0s=[[10.0]] * J),
)

# Extract marginal posteriors
hyper_result  = hproblem.extract_hyper(result)
group_results = [hproblem.extract_group(result, j) for j in range(J)]

HierarchicalProblem

Joint Bayesian inference over hyperparameters and group-level parameters.

Models the situation where you have J groups (e.g. J nominally identical structures) each with local parameters :math:\theta_j, linked through shared hyperparameters :math:\phi:

.. math::

p(\phi) \quad\text{(hyperprior)}

p(\theta_j \mid \phi) \quad\text{(group-level prior)}

p(y_j \mid \theta_j) \quad\text{(per-group likelihood)}

The joint log-posterior is:

.. math::

\log p(\phi, \theta_1,\ldots,\theta_J \mid y)
  = \log p(\phi)
  + \sum_j \log p(\theta_j \mid \phi)
  + \sum_j \log p(y_j \mid \theta_j)

The full parameter vector is laid out as:

.. code-block::

theta = [phi_0, ..., phi_{n_hyper-1},
         theta_1_0, ..., theta_1_{n_group-1},
         ...
         theta_J_0, ..., theta_J_{n_group-1}]

This object exposes the same interface as :class:Problem and is therefore compatible with all samplers.

Parameters:

Name Type Description Default
hyperprior callable

log p(phi) -> float. phi has shape (n_hyper,).

 required
group_prior callable

log p(theta_j | phi) -> float. Called as group_prior(theta_j, phi).

 required
group_likelihoods list of callable

One callable per group: log p(y_j | theta_j) -> float. Must be a list, one entry per group; the number of groups is taken from its length. A bare callable is rejected rather than broadcast, since sharing one likelihood across every group is almost always a mistake.

 required
n_hyper int

Dimension of phi.

 required
n_group int

Dimension of each group's local parameter vector theta_j.

 required
param_names_hyper list of str

Names for the hyperparameters.

 None
param_names_group list of str

Names for one group's parameters; automatically suffixed with the group index, e.g. ["k[0]", "k[1]", ...].

 None

Examples:

5 structures sharing a stiffness population::

def hyperprior(phi):
    mu, log_s = phi
    return -0.5 * ((mu - 10) / 3)**2 - 0.5 * log_s**2

def group_prior(theta_j, phi):
    mu, log_s = phi
    sigma = np.exp(log_s)
    return -0.5 * ((theta_j[0] - mu) / sigma)**2 - log_s

likelihoods = [make_likelihood(y_j) for y_j in datasets]

hproblem = mc.HierarchicalProblem(
    hyperprior=hyperprior,
    group_prior=group_prior,
    group_likelihoods=likelihoods,
    n_hyper=2,
    n_group=1,
    param_names_hyper=["mu_k", "log_sigma_k"],
    param_names_group=["k"],
)

result = mc.DRAM(n_samples=20_000, initial_cov=np.eye(hproblem.n_params)).run(
    hproblem, x0=hproblem.default_x0([10.0, np.log(1.0)],
                                      [[10.0]] * 5))

hyper_result  = hproblem.extract_hyper(result)
group_results = [hproblem.extract_group(result, j) for j in range(5)]
Source code in mcmckit/core/hierarchical.py 
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
class HierarchicalProblem:
    """Joint Bayesian inference over hyperparameters and group-level parameters.

    Models the situation where you have J groups (e.g. J nominally identical
    structures) each with local parameters :math:`\\theta_j`, linked through
    shared hyperparameters :math:`\\phi`:

    .. math::

        p(\\phi) \\quad\\text{(hyperprior)}

        p(\\theta_j \\mid \\phi) \\quad\\text{(group-level prior)}

        p(y_j \\mid \\theta_j) \\quad\\text{(per-group likelihood)}

    The joint log-posterior is:

    .. math::

        \\log p(\\phi, \\theta_1,\\ldots,\\theta_J \\mid y)
          = \\log p(\\phi)
          + \\sum_j \\log p(\\theta_j \\mid \\phi)
          + \\sum_j \\log p(y_j \\mid \\theta_j)

    The full parameter vector is laid out as:

    .. code-block::

        theta = [phi_0, ..., phi_{n_hyper-1},
                 theta_1_0, ..., theta_1_{n_group-1},
                 ...
                 theta_J_0, ..., theta_J_{n_group-1}]

    This object exposes the same interface as :class:`Problem` and is
    therefore compatible with all samplers.

    Parameters
    ----------
    hyperprior : callable
        ``log p(phi) -> float``.  ``phi`` has shape ``(n_hyper,)``.
    group_prior : callable
        ``log p(theta_j | phi) -> float``.  Called as
        ``group_prior(theta_j, phi)``.
    group_likelihoods : list of callable
        One callable per group: ``log p(y_j | theta_j) -> float``.
        Must be a list, one entry per group; the number of groups is taken from
        its length. A bare callable is rejected rather than broadcast, since
        sharing one likelihood across every group is almost always a mistake.
    n_hyper : int
        Dimension of ``phi``.
    n_group : int
        Dimension of each group's local parameter vector ``theta_j``.
    param_names_hyper : list of str, optional
        Names for the hyperparameters.
    param_names_group : list of str, optional
        Names for one group's parameters; automatically suffixed with the
        group index, e.g. ``["k[0]", "k[1]", ...]``.

    Examples
    --------
    5 structures sharing a stiffness population::

        def hyperprior(phi):
            mu, log_s = phi
            return -0.5 * ((mu - 10) / 3)**2 - 0.5 * log_s**2

        def group_prior(theta_j, phi):
            mu, log_s = phi
            sigma = np.exp(log_s)
            return -0.5 * ((theta_j[0] - mu) / sigma)**2 - log_s

        likelihoods = [make_likelihood(y_j) for y_j in datasets]

        hproblem = mc.HierarchicalProblem(
            hyperprior=hyperprior,
            group_prior=group_prior,
            group_likelihoods=likelihoods,
            n_hyper=2,
            n_group=1,
            param_names_hyper=["mu_k", "log_sigma_k"],
            param_names_group=["k"],
        )

        result = mc.DRAM(n_samples=20_000, initial_cov=np.eye(hproblem.n_params)).run(
            hproblem, x0=hproblem.default_x0([10.0, np.log(1.0)],
                                              [[10.0]] * 5))

        hyper_result  = hproblem.extract_hyper(result)
        group_results = [hproblem.extract_group(result, j) for j in range(5)]
    """

    def __init__(self, hyperprior, group_prior, group_likelihoods,
                 n_hyper, n_group,
                 param_names_hyper=None, param_names_group=None):

        self._hyperprior = hyperprior
        self._group_prior = group_prior

        # Accept a single likelihood shared by all groups
        if callable(group_likelihoods) and not isinstance(group_likelihoods, list):
            raise ValueError(
                "group_likelihoods must be a list of callables, one per group."
            )
        self._group_likelihoods = list(group_likelihoods)
        self.n_groups = len(self._group_likelihoods)
        self.n_hyper = n_hyper
        self.n_group = n_group

        # Build param_names
        hyper_names = (list(param_names_hyper)
                       if param_names_hyper is not None
                       else [f"phi[{i}]" for i in range(n_hyper)])
        if param_names_group is not None:
            base = list(param_names_group)
        else:
            base = [f"theta[{i}]" for i in range(n_group)]

        group_names = [f"{name}_{j}"
                       for j in range(self.n_groups)
                       for name in base]
        self.param_names = hyper_names + group_names

    # ------------------------------------------------------------------
    # Helpers: split / assemble the full parameter vector
    # ------------------------------------------------------------------

    @property
    def n_params(self):
        """Total number of parameters: n_hyper + n_groups * n_group."""
        return self.n_hyper + self.n_groups * self.n_group

    def split(self, theta):
        """Split full theta into (phi, [theta_0, ..., theta_{J-1}]).

        Parameters
        ----------
        theta : array-like, shape (n_params,)

        Returns
        -------
        phi : np.ndarray, shape (n_hyper,)
        group_thetas : list of np.ndarray, each shape (n_group,)
        """
        theta = np.asarray(theta, dtype=float)
        phi = theta[:self.n_hyper]
        group_thetas = [
            theta[self.n_hyper + j * self.n_group:
                  self.n_hyper + (j + 1) * self.n_group]
            for j in range(self.n_groups)
        ]
        return phi, group_thetas

    def default_x0(self, phi0, group_x0s):
        """Assemble a starting point from hyper and group initial values.

        Parameters
        ----------
        phi0 : array-like, shape (n_hyper,)
        group_x0s : list of array-like, each shape (n_group,)
            One entry per group.

        Returns
        -------
        np.ndarray, shape (n_params,)
        """
        parts = [np.asarray(phi0, dtype=float)]
        for x in group_x0s:
            parts.append(np.asarray(x, dtype=float))
        return np.concatenate(parts)

    # ------------------------------------------------------------------
    # Problem interface (compatible with all samplers)
    # ------------------------------------------------------------------

    def log_prior(self, theta):
        """log p(phi) + sum_j log p(theta_j | phi)."""
        phi, group_thetas = self.split(theta)
        lp = float(self._hyperprior(phi))
        if not np.isfinite(lp):
            return -np.inf
        for theta_j in group_thetas:
            lp += float(self._group_prior(theta_j, phi))
            if not np.isfinite(lp):
                return -np.inf
        return lp

    def log_likelihood(self, theta):
        """sum_j log p(y_j | theta_j)."""
        _, group_thetas = self.split(theta)
        ll = 0.0
        for theta_j, lik_j in zip(group_thetas, self._group_likelihoods):
            ll += float(lik_j(theta_j))
            if not np.isfinite(ll):
                return -np.inf
        return ll

    def log_posterior(self, theta):
        """log p(phi, theta_1,...,theta_J | y)."""
        lp = self.log_prior(theta)
        if not np.isfinite(lp):
            return -np.inf
        return lp + self.log_likelihood(theta)

    # ------------------------------------------------------------------
    # Extract marginal Results after sampling
    # ------------------------------------------------------------------

    def extract_hyper(self, result):
        """Return a Result containing only the hyperparameter samples.

        Parameters
        ----------
        result : Result
            Full posterior result from running a sampler on this problem.

        Returns
        -------
        Result, shape (n_samples, n_hyper)
        """
        names = self.param_names[:self.n_hyper]
        return Result(
            samples=result.samples[:, :self.n_hyper],
            log_posteriors=result.log_posteriors,
            param_names=names,
        )

    def extract_group(self, result, j):
        """Return a Result containing only group j's parameter samples.

        Parameters
        ----------
        result : Result
        j : int
            Group index (0-based).

        Returns
        -------
        Result, shape (n_samples, n_group)
        """
        if j < 0 or j >= self.n_groups:
            raise ValueError(f"j must be in [0, {self.n_groups - 1}], got {j}.")
        start = self.n_hyper + j * self.n_group
        end   = start + self.n_group
        # recover group param names from full param_names
        names = self.param_names[start:end]
        return Result(
            samples=result.samples[:, start:end],
            log_posteriors=result.log_posteriors,
            param_names=names,
        )

Attributes

n_params property 

n_params

Total number of parameters: n_hyper + n_groups * n_group.

Methods:

split 

split(theta)

Split full theta into (phi, [theta_0, ..., theta_{J-1}]).

Parameters:

Name Type Description Default
theta (array - like, shape(n_params))
required

Returns:

Name Type Description
phi (ndarray, shape(n_hyper))
group_thetas list of np.ndarray, each shape (n_group,)
Source code in mcmckit/core/hierarchical.py 
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
def split(self, theta):
    """Split full theta into (phi, [theta_0, ..., theta_{J-1}]).

    Parameters
    ----------
    theta : array-like, shape (n_params,)

    Returns
    -------
    phi : np.ndarray, shape (n_hyper,)
    group_thetas : list of np.ndarray, each shape (n_group,)
    """
    theta = np.asarray(theta, dtype=float)
    phi = theta[:self.n_hyper]
    group_thetas = [
        theta[self.n_hyper + j * self.n_group:
              self.n_hyper + (j + 1) * self.n_group]
        for j in range(self.n_groups)
    ]
    return phi, group_thetas

default_x0 

default_x0(phi0, group_x0s)

Assemble a starting point from hyper and group initial values.

Parameters:

Name Type Description Default
phi0 (array - like, shape(n_hyper))
required
group_x0s list of array-like, each shape (n_group,)

One entry per group.

 required

Returns:

Type Description
(ndarray, shape(n_params))
Source code in mcmckit/core/hierarchical.py 
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
def default_x0(self, phi0, group_x0s):
    """Assemble a starting point from hyper and group initial values.

    Parameters
    ----------
    phi0 : array-like, shape (n_hyper,)
    group_x0s : list of array-like, each shape (n_group,)
        One entry per group.

    Returns
    -------
    np.ndarray, shape (n_params,)
    """
    parts = [np.asarray(phi0, dtype=float)]
    for x in group_x0s:
        parts.append(np.asarray(x, dtype=float))
    return np.concatenate(parts)

log_prior 

log_prior(theta)

log p(phi) + sum_j log p(theta_j | phi).

Source code in mcmckit/core/hierarchical.py 
181
182
183
184
185
186
187
188
189
190
191
def log_prior(self, theta):
    """log p(phi) + sum_j log p(theta_j | phi)."""
    phi, group_thetas = self.split(theta)
    lp = float(self._hyperprior(phi))
    if not np.isfinite(lp):
        return -np.inf
    for theta_j in group_thetas:
        lp += float(self._group_prior(theta_j, phi))
        if not np.isfinite(lp):
            return -np.inf
    return lp

log_likelihood 

log_likelihood(theta)

sum_j log p(y_j | theta_j).

Source code in mcmckit/core/hierarchical.py 
193
194
195
196
197
198
199
200
201
def log_likelihood(self, theta):
    """sum_j log p(y_j | theta_j)."""
    _, group_thetas = self.split(theta)
    ll = 0.0
    for theta_j, lik_j in zip(group_thetas, self._group_likelihoods):
        ll += float(lik_j(theta_j))
        if not np.isfinite(ll):
            return -np.inf
    return ll

log_posterior 

log_posterior(theta)

log p(phi, theta_1,...,theta_J | y).

Source code in mcmckit/core/hierarchical.py 
203
204
205
206
207
208
def log_posterior(self, theta):
    """log p(phi, theta_1,...,theta_J | y)."""
    lp = self.log_prior(theta)
    if not np.isfinite(lp):
        return -np.inf
    return lp + self.log_likelihood(theta)

extract_hyper 

extract_hyper(result)

Return a Result containing only the hyperparameter samples.

Parameters:

Name Type Description Default
result Result

Full posterior result from running a sampler on this problem.

 required

Returns:

Type Description
(Result, shape(n_samples, n_hyper))
Source code in mcmckit/core/hierarchical.py 
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
def extract_hyper(self, result):
    """Return a Result containing only the hyperparameter samples.

    Parameters
    ----------
    result : Result
        Full posterior result from running a sampler on this problem.

    Returns
    -------
    Result, shape (n_samples, n_hyper)
    """
    names = self.param_names[:self.n_hyper]
    return Result(
        samples=result.samples[:, :self.n_hyper],
        log_posteriors=result.log_posteriors,
        param_names=names,
    )

extract_group 

extract_group(result, j)

Return a Result containing only group j's parameter samples.

Parameters:

Name Type Description Default
result Result
required
j int

Group index (0-based).

 required

Returns:

Type Description
(Result, shape(n_samples, n_group))
Source code in mcmckit/core/hierarchical.py 
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
def extract_group(self, result, j):
    """Return a Result containing only group j's parameter samples.

    Parameters
    ----------
    result : Result
    j : int
        Group index (0-based).

    Returns
    -------
    Result, shape (n_samples, n_group)
    """
    if j < 0 or j >= self.n_groups:
        raise ValueError(f"j must be in [0, {self.n_groups - 1}], got {j}.")
    start = self.n_hyper + j * self.n_group
    end   = start + self.n_group
    # recover group param names from full param_names
    names = self.param_names[start:end]
    return Result(
        samples=result.samples[:, start:end],
        log_posteriors=result.log_posteriors,
        param_names=names,
    )