Getting started with parameteric priors¶
Here we introduce how to specify the elicitation method for a parametric prior.
Imports¶
In [1]:
Copied!
import os
os.environ["TF_ENABLE_ONEDNN_OPTS"] = "0"
import copy
from typing import Any
import numpy as np
import tensorflow as tf
import tensorflow_probability as tfp
import elicito as el
tfd = tfp.distributions
import os
os.environ["TF_ENABLE_ONEDNN_OPTS"] = "0"
import copy
from typing import Any
import numpy as np
import tensorflow as tf
import tensorflow_probability as tfp
import elicito as el
tfd = tfp.distributions
2025-11-12 14:29:45.530206: I external/local_xla/xla/tsl/cuda/cudart_stub.cc:31] Could not find cuda drivers on your machine, GPU will not be used. 2025-11-12 14:29:45.577460: I tensorflow/core/platform/cpu_feature_guard.cc:210] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations. To enable the following instructions: AVX2 FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags.
2025-11-12 14:29:47.255140: I external/local_xla/xla/tsl/cuda/cudart_stub.cc:31] Could not find cuda drivers on your machine, GPU will not be used.
2025-11-12 14:29:48.940254: E external/local_xla/xla/stream_executor/cuda/cuda_platform.cc:51] failed call to cuInit: INTERNAL: CUDA error: Failed call to cuInit: UNKNOWN ERROR (303)
The Model¶
Probabilistic model¶
$$ \begin{align*} \beta_0 &\sim \text{Normal}(\mu_0, \sigma_0) \\ \beta_1 &\sim \text{Normal}(\mu_1, \sigma_1) \\ \sigma &\sim \text{HalfNormal}(\sigma_2) \\ \mu &= \beta_0 + \beta_1X \\ y_{pred} &\sim \text{Normal}(\mu, \sigma) \end{align*} $$
Implementation¶
Predictor¶
In [2]:
Copied!
# create a predictor ranging from 1 to 200
# standardize predictor
# select the 25th (X0), 50th (X1), and 75th (X2) quantile of the std.
# predictor for querying the expert
def X_design(N: int, quantiles: list[float]) -> np.ndarray:
"""
Compute design matrix
Parameters
----------
N
number of observations
quantiles
list of quantiles
Returns
-------
:
design matrix
"""
X = tf.cast(np.arange(N), tf.float32)
X_std = X / tf.math.reduce_std(X)
X_sel = tfp.stats.percentile(X_std, quantiles)
X_design = tf.stack([tf.ones(X_sel.shape), X_sel], -1)
return X_design
X_design(N=30, quantiles=[25, 50, 75])
# create a predictor ranging from 1 to 200
# standardize predictor
# select the 25th (X0), 50th (X1), and 75th (X2) quantile of the std.
# predictor for querying the expert
def X_design(N: int, quantiles: list[float]) -> np.ndarray:
"""
Compute design matrix
Parameters
----------
N
number of observations
quantiles
list of quantiles
Returns
-------
:
design matrix
"""
X = tf.cast(np.arange(N), tf.float32)
X_std = X / tf.math.reduce_std(X)
X_sel = tfp.stats.percentile(X_std, quantiles)
X_design = tf.stack([tf.ones(X_sel.shape), X_sel], -1)
return X_design
X_design(N=30, quantiles=[25, 50, 75])
Out[2]:
<tf.Tensor: shape=(3, 2), dtype=float32, numpy=
array([[1. , 0.80873984],
[1. , 1.6174797 ],
[1. , 2.5417538 ]], dtype=float32)>
Generative model¶
In [3]:
Copied!
class ToyModel:
"""
generative model
"""
def __call__(self, prior_samples: Any, design_matrix: Any) -> dict[str, Any]:
"""
Compute target quantities from generative model
Parameters
----------
prior_samples
prior samples
design_matrix
design matrix
Returns
-------
:
dictionary with target quantities
"""
# linear predictor
epred = tf.matmul(prior_samples[:, :, :-1], design_matrix, transpose_b=True)
# data-generating model
likelihood = tfd.Normal(
loc=epred, scale=tf.expand_dims(prior_samples[:, :, -1], -1)
)
# prior predictive distribution
ypred = likelihood.sample()
# selected observations
y_X0, y_X1, y_X2 = (ypred[:, :, 0], ypred[:, :, 1], ypred[:, :, 2])
return dict(y_X0=y_X0, y_X1=y_X1, y_X2=y_X2, ypred=ypred, epred=epred)
class ToyModel:
"""
generative model
"""
def __call__(self, prior_samples: Any, design_matrix: Any) -> dict[str, Any]:
"""
Compute target quantities from generative model
Parameters
----------
prior_samples
prior samples
design_matrix
design matrix
Returns
-------
:
dictionary with target quantities
"""
# linear predictor
epred = tf.matmul(prior_samples[:, :, :-1], design_matrix, transpose_b=True)
# data-generating model
likelihood = tfd.Normal(
loc=epred, scale=tf.expand_dims(prior_samples[:, :, -1], -1)
)
# prior predictive distribution
ypred = likelihood.sample()
# selected observations
y_X0, y_X1, y_X2 = (ypred[:, :, 0], ypred[:, :, 1], ypred[:, :, 2])
return dict(y_X0=y_X0, y_X1=y_X1, y_X2=y_X2, ypred=ypred, epred=epred)
Generative model¶
In [4]:
Copied!
# specify the model
model = el.model(obj=ToyModel, design_matrix=X_design(N=30, quantiles=[25, 50, 75]))
# specify the model
model = el.model(obj=ToyModel, design_matrix=X_design(N=30, quantiles=[25, 50, 75]))
Model parameters¶
- intercept with normal prior $\beta_0$
- slope with normal prior $\beta_1$
- random noise with halfnormal prior $\sigma$
To be learned hyperparameters $ \lambda = (\mu_0, \sigma_0, \mu_1, \sigma_1, \sigma_2) $
- scale parameters ($\sigma_0, \sigma_1, \sigma_2$) are constrained to be positive
In [5]:
Copied!
parameters = [
el.parameter(
name="beta0",
family=tfd.Normal,
hyperparams=dict(loc=el.hyper("mu0"), scale=el.hyper("sigma0", lower=0)),
),
el.parameter(
name="beta1",
family=tfd.Normal,
hyperparams=dict(loc=el.hyper("mu1"), scale=el.hyper("sigma1", lower=0)),
),
el.parameter(
name="sigma",
family=tfd.HalfNormal,
hyperparams=dict(scale=el.hyper("sigma2", lower=0)),
),
]
parameters = [
el.parameter(
name="beta0",
family=tfd.Normal,
hyperparams=dict(loc=el.hyper("mu0"), scale=el.hyper("sigma0", lower=0)),
),
el.parameter(
name="beta1",
family=tfd.Normal,
hyperparams=dict(loc=el.hyper("mu1"), scale=el.hyper("sigma1", lower=0)),
),
el.parameter(
name="sigma",
family=tfd.HalfNormal,
hyperparams=dict(scale=el.hyper("sigma2", lower=0)),
),
]
Target quantities and elicitation techniques¶
Target quantities
- query expert regarding prior predictions $y \mid X_{i}$ with $i$ being the 25th, 50th, and 75th quantile of the predictor.
Elicitation technique
- query each observation using quantile-based elicitation using $Q_p(y \mid X)$ for $p=25, 50, 75$
Specifying discrepancy measure and weight for single loss components
- prior predictions $y \mid X_i$:
- discrepancy measure: Maximum Mean Discrepancy with energy kernel
- weight of loss component: 1.0
- $R^2$:
- L2 loss
- weight of loss component: 0.5
In [6]:
Copied!
def custom_r2(ypred, epred):
"""Compute coefficient of determination"""
var_epred = tf.math.reduce_variance(epred, -1)
# variance of difference between ypred and epred
var_diff = tf.math.reduce_variance(tf.subtract(ypred, epred), -1)
var_total = var_epred + var_diff
# variance of linear predictor divided by total variance
return tf.divide(var_epred, var_total)
targets = [
el.target(
name="y_X0",
query=el.queries.quantiles((0.05, 0.25, 0.50, 0.75, 0.95)),
loss=el.losses.MMD2(kernel="energy"),
weight=1.0,
),
el.target(
name="y_X1",
query=el.queries.quantiles((0.05, 0.25, 0.50, 0.75, 0.95)),
loss=el.losses.MMD2(kernel="energy"),
weight=1.0,
),
el.target(
name="y_X2",
query=el.queries.quantiles((0.05, 0.25, 0.50, 0.75, 0.95)),
loss=el.losses.MMD2(kernel="energy"),
weight=1.0,
),
el.target(
name="R2",
query=el.queries.quantiles((0.05, 0.50, 0.95)),
loss=el.losses.L2,
weight=0.5,
target_method=custom_r2,
),
]
def custom_r2(ypred, epred):
"""Compute coefficient of determination"""
var_epred = tf.math.reduce_variance(epred, -1)
# variance of difference between ypred and epred
var_diff = tf.math.reduce_variance(tf.subtract(ypred, epred), -1)
var_total = var_epred + var_diff
# variance of linear predictor divided by total variance
return tf.divide(var_epred, var_total)
targets = [
el.target(
name="y_X0",
query=el.queries.quantiles((0.05, 0.25, 0.50, 0.75, 0.95)),
loss=el.losses.MMD2(kernel="energy"),
weight=1.0,
),
el.target(
name="y_X1",
query=el.queries.quantiles((0.05, 0.25, 0.50, 0.75, 0.95)),
loss=el.losses.MMD2(kernel="energy"),
weight=1.0,
),
el.target(
name="y_X2",
query=el.queries.quantiles((0.05, 0.25, 0.50, 0.75, 0.95)),
loss=el.losses.MMD2(kernel="energy"),
weight=1.0,
),
el.target(
name="R2",
query=el.queries.quantiles((0.05, 0.50, 0.95)),
loss=el.losses.L2,
weight=0.5,
target_method=custom_r2,
),
]
Expert elicitation¶
Instead of querying a "real" expert, we define a ground truth (i.e., oracle) and simulate the oracle-elicited statistics
In [7]:
Copied!
# specify ground truth
ground_truth = {
"beta0": tfd.Normal(loc=5, scale=1),
"beta1": tfd.Normal(loc=2, scale=1),
"sigma": tfd.HalfNormal(scale=7.0),
}
# define oracle
expert = el.expert.simulator(ground_truth=ground_truth, num_samples=10_000)
# specify ground truth
ground_truth = {
"beta0": tfd.Normal(loc=5, scale=1),
"beta1": tfd.Normal(loc=2, scale=1),
"sigma": tfd.HalfNormal(scale=7.0),
}
# define oracle
expert = el.expert.simulator(ground_truth=ground_truth, num_samples=10_000)
Training¶
Learn prior distributions based on expert data
In [8]:
Copied!
eliobj = el.Elicit(
model=model,
parameters=parameters,
targets=targets,
expert=expert,
optimizer=el.optimizer(
optimizer=tf.keras.optimizers.Adam, learning_rate=0.1, clipnorm=1.0
),
trainer=el.trainer(method="parametric_prior", seed=2025, epochs=500, progress=0),
initializer=el.initializer(
method="sobol",
iterations=32,
distribution=el.initialization.uniform(radius=2.0, mean=0.0),
),
)
eliobj = el.Elicit(
model=model,
parameters=parameters,
targets=targets,
expert=expert,
optimizer=el.optimizer(
optimizer=tf.keras.optimizers.Adam, learning_rate=0.1, clipnorm=1.0
),
trainer=el.trainer(method="parametric_prior", seed=2025, epochs=500, progress=0),
initializer=el.initializer(
method="sobol",
iterations=32,
distribution=el.initialization.uniform(radius=2.0, mean=0.0),
),
)
Inspect summary of eliobj
In [9]:
Copied!
eliobj
eliobj
Out[9]:
Model hyperparameters: 5 Model parameters: 3 Targets -> Elicited summaries (loss components): 4 - y_X0 (128, 200) -> quantiles_y_X0 (128, 5) - y_X1 (128, 200) -> quantiles_y_X1 (128, 5) - y_X2 (128, 200) -> quantiles_y_X2 (128, 5) - R2 (128, 200) -> quantiles_R2 (128, 3) Prior samples: 200 (128, 200, 3) Batch size: 128 Epochs: 500 Method: parametric_prior Seed: 2025 Optimizer: Adam(lr=0.1) Initializer: (method: sobol, iterations: 32)
Fit eliobj
In [10]:
Copied!
eliobj.fit(parallel=el.utils.parallel())
eliobj.fit(parallel=el.utils.parallel())
2025-11-12 14:29:50.664487: I external/local_xla/xla/tsl/cuda/cudart_stub.cc:31] Could not find cuda drivers on your machine, GPU will not be used. 2025-11-12 14:29:50.757890: I external/local_xla/xla/tsl/cuda/cudart_stub.cc:31] Could not find cuda drivers on your machine, GPU will not be used. 2025-11-12 14:29:50.783254: I external/local_xla/xla/tsl/cuda/cudart_stub.cc:31] Could not find cuda drivers on your machine, GPU will not be used. 2025-11-12 14:29:50.794947: I tensorflow/core/platform/cpu_feature_guard.cc:210] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations. To enable the following instructions: AVX2 FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags. 2025-11-12 14:29:50.827570: I external/local_xla/xla/tsl/cuda/cudart_stub.cc:31] Could not find cuda drivers on your machine, GPU will not be used.
2025-11-12 14:29:50.906301: I tensorflow/core/platform/cpu_feature_guard.cc:210] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations. To enable the following instructions: AVX2 FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags. 2025-11-12 14:29:50.933578: I tensorflow/core/platform/cpu_feature_guard.cc:210] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations. To enable the following instructions: AVX2 FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags. 2025-11-12 14:29:50.979509: I tensorflow/core/platform/cpu_feature_guard.cc:210] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations. To enable the following instructions: AVX2 FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags.
2025-11-12 14:29:56.072558: I external/local_xla/xla/tsl/cuda/cudart_stub.cc:31] Could not find cuda drivers on your machine, GPU will not be used. 2025-11-12 14:29:56.248407: I external/local_xla/xla/tsl/cuda/cudart_stub.cc:31] Could not find cuda drivers on your machine, GPU will not be used.
2025-11-12 14:29:56.399707: I external/local_xla/xla/tsl/cuda/cudart_stub.cc:31] Could not find cuda drivers on your machine, GPU will not be used. 2025-11-12 14:29:56.404394: I external/local_xla/xla/tsl/cuda/cudart_stub.cc:31] Could not find cuda drivers on your machine, GPU will not be used.
2025-11-12 14:30:02.051841: E external/local_xla/xla/stream_executor/cuda/cuda_platform.cc:51] failed call to cuInit: INTERNAL: CUDA error: Failed call to cuInit: UNKNOWN ERROR (303)
2025-11-12 14:30:02.384746: E external/local_xla/xla/stream_executor/cuda/cuda_platform.cc:51] failed call to cuInit: INTERNAL: CUDA error: Failed call to cuInit: UNKNOWN ERROR (303)
2025-11-12 14:30:02.595440: E external/local_xla/xla/stream_executor/cuda/cuda_platform.cc:51] failed call to cuInit: INTERNAL: CUDA error: Failed call to cuInit: UNKNOWN ERROR (303) 2025-11-12 14:30:02.633780: E external/local_xla/xla/stream_executor/cuda/cuda_platform.cc:51] failed call to cuInit: INTERNAL: CUDA error: Failed call to cuInit: UNKNOWN ERROR (303)
Inspect results
In [11]:
Copied!
eliobj.results
eliobj.results
Out[11]:
<xarray.DatasetView> Size: 0B
Dimensions: ()
Data variables:
*empty*In [12]:
Copied!
el.plots.initialization(eliobj,
titles=["$\mu_0$", "$\sigma_0$", "$\mu_2$", "$\sigma_1$", "$\sigma_2$"],
cols=4);
el.plots.initialization(eliobj,
titles=["$\mu_0$", "$\sigma_0$", "$\mu_2$", "$\sigma_1$", "$\sigma_2$"],
cols=4);
Convergence - Loss¶
In [13]:
Copied!
el.plots.loss(eliobj);
el.plots.loss(eliobj);
Convergence - hyperparameters¶
In [14]:
Copied!
el.plots.hyperparameter(eliobj,
titles=["$\mu_0$", "$\sigma_0$", "$\mu_2$", "$\sigma_1$", "$\sigma_2$"],
cols=5);
el.plots.hyperparameter(eliobj,
titles=["$\mu_0$", "$\sigma_0$", "$\mu_2$", "$\sigma_1$", "$\sigma_2$"],
cols=5);
Expert expectations¶
In [15]:
Copied!
el.plots.elicits(eliobj, cols=4);
el.plots.elicits(eliobj, cols=4);
Learned prior distributions¶
In [16]:
Copied!
el.plots.prior_marginals(eliobj, cols=3);
el.plots.prior_marginals(eliobj, cols=3);
In [17]:
Copied!
el.plots.priorpredictive(eliobj, target="R2");
el.plots.priorpredictive(eliobj, target="R2");
In [18]:
Copied!
el.plots.prior_averaging(eliobj);
el.plots.prior_averaging(eliobj);
INFO: Reset cols=3
INFO: Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting.
INFO: Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting.
INFO: Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting.
Add-on: Shared parameters¶
In [19]:
Copied!
# create a copy of eliobj
eliobj_shared = copy.deepcopy(eliobj)
# share sigma hyperparameter of intercept and slope parameter
parameters_shared = [
el.parameter(
name="beta0",
family=tfd.Normal,
hyperparams=dict(
loc=el.hyper("mu0"), scale=el.hyper("sigma1", lower=0, shared=True)
),
),
el.parameter(
name="beta1",
family=tfd.Normal,
hyperparams=dict(
loc=el.hyper("mu1"), scale=el.hyper("sigma1", lower=0, shared=True)
),
),
el.parameter(
name="sigma",
family=tfd.HalfNormal,
hyperparams=dict(scale=el.hyper("sigma2", lower=0)),
),
]
# update parameters in eliobj
eliobj_shared.update(parameters=parameters_shared)
# refit the model
eliobj_shared.fit()
# create a copy of eliobj
eliobj_shared = copy.deepcopy(eliobj)
# share sigma hyperparameter of intercept and slope parameter
parameters_shared = [
el.parameter(
name="beta0",
family=tfd.Normal,
hyperparams=dict(
loc=el.hyper("mu0"), scale=el.hyper("sigma1", lower=0, shared=True)
),
),
el.parameter(
name="beta1",
family=tfd.Normal,
hyperparams=dict(
loc=el.hyper("mu1"), scale=el.hyper("sigma1", lower=0, shared=True)
),
),
el.parameter(
name="sigma",
family=tfd.HalfNormal,
hyperparams=dict(scale=el.hyper("sigma2", lower=0)),
),
]
# update parameters in eliobj
eliobj_shared.update(parameters=parameters_shared)
# refit the model
eliobj_shared.fit()
INFO: Results have been reset.
In [20]:
Copied!
el.plots.initialization(eliobj_shared, cols=4);
el.plots.initialization(eliobj_shared, cols=4);
Convergence - Loss¶
In [21]:
Copied!
el.plots.loss(eliobj_shared);
el.plots.loss(eliobj_shared);
Convergence - Hyperparameters¶
In [22]:
Copied!
el.plots.hyperparameter(eliobj_shared);
el.plots.hyperparameter(eliobj_shared);
Elicited statistics¶
In [23]:
Copied!
el.plots.elicits(eliobj_shared, cols=4);
el.plots.elicits(eliobj_shared, cols=4);
Learned prior distributions¶
In [24]:
Copied!
el.plots.prior_marginals(eliobj_shared);
el.plots.prior_marginals(eliobj_shared);
INFO: Reset cols=3
Add-on: Use expert data as input¶
In [25]:
Copied!
# create a copy of eliobj
eliobj_dat = copy.deepcopy(eliobj)
# use dictionary of elicited expert data (instead of simulating data)
expert_dat = {
"quantiles_y_X0": [-12.5, -0.6, 3.3, 7.1, 19.1],
"quantiles_y_X1": [-11.2, 1.5, 5.0, 8.8, 20.4],
"quantiles_y_X2": [-9.3, 3.1, 6.8, 10.5, 23.3],
"quantiles_R2": [0.001, 0.09, 0.96]
}
# update expert in eliobj
eliobj_dat.update(expert=el.expert.data(dat=expert_dat))
# refit the model
eliobj_dat.fit()
# create a copy of eliobj
eliobj_dat = copy.deepcopy(eliobj)
# use dictionary of elicited expert data (instead of simulating data)
expert_dat = {
"quantiles_y_X0": [-12.5, -0.6, 3.3, 7.1, 19.1],
"quantiles_y_X1": [-11.2, 1.5, 5.0, 8.8, 20.4],
"quantiles_y_X2": [-9.3, 3.1, 6.8, 10.5, 23.3],
"quantiles_R2": [0.001, 0.09, 0.96]
}
# update expert in eliobj
eliobj_dat.update(expert=el.expert.data(dat=expert_dat))
# refit the model
eliobj_dat.fit()
INFO: Results have been reset.
In [26]:
Copied!
el.plots.initialization(eliobj_dat, cols=5);
el.plots.initialization(eliobj_dat, cols=5);
Convergence - Loss¶
In [27]:
Copied!
el.plots.loss(eliobj_dat);
el.plots.loss(eliobj_dat);
Convergence - Hyperparameters¶
In [28]:
Copied!
el.plots.hyperparameter(eliobj_dat, cols=5);
el.plots.hyperparameter(eliobj_dat, cols=5);
Elicited statistics¶
In [29]:
Copied!
el.plots.elicits(eliobj_dat, cols=4);
el.plots.elicits(eliobj_dat, cols=4);
Learned prior distributions¶
In [30]:
Copied!
el.plots.prior_marginals(eliobj_dat);
el.plots.prior_marginals(eliobj_dat);
INFO: Reset cols=3
