Adds support for multiple targets with weights to ppe

preliz/internal/optimization.py

@@ -227,12 +227,17 @@ def optimize_beta_mode(lower, upper, tau_not, mode, dist, mass, prob):
             tau_not += 0.5 * tau_not
 
 
-def optimize_pymc_model(fmodel, target, draws, prior, initial_guess, bounds, var_info, p_model):
+def optimize_pymc_model(
+    fmodel, target, draws, prior, initial_guess, bounds, var_info, p_model, rng
+):
     for _ in range(400):
         # can we sample systematically from these and less random?
         # This should be more flexible and allow other targets than just
         # a preliz distribution
-        obs = target.rvs(draws, random_state=np.random.default_rng(_))
+        if isinstance(target, list):
+            obs = get_weighted_rvs(target, draws, rng)
+        else:
+            obs = target.rvs(draws, random_state=rng)
         result = minimize(
             fmodel,
             initial_guess,
@@ -452,3 +457,11 @@ def func(x, dist, q):
         left, right = right, right * factor
 
     return brentq(func, left, right, args=(dist, q))
+
+
+def get_weighted_rvs(target, size, rng):
+    targets = [t[0] for t in target]
+    weights = [t[1] for t in target]
+    target_rnd_choices = np.random.choice(len(targets), size=size, p=weights)
+    samples = [target.rvs(size, random_state=rng) for target in targets]
+    return np.choose(target_rnd_choices, samples)

preliz/predictive/ppe.py

@@ -1,6 +1,7 @@
 """Projective predictive elicitation."""
 
 import logging
+import numpy as np
 
 from preliz.internal.optimization import optimize_pymc_model
 from preliz.ppls.pymc_io import (
@@ -15,7 +16,7 @@
 _log = logging.getLogger("preliz")
 
 
-def ppe(model, target):
+def ppe(model, target, seed=0):
     """
     Projective Predictive Elicitation.
 
@@ -26,13 +27,15 @@ def ppe(model, target):
     ----------
     model : a probabilistic model
         Currently it only works with PyMC model. More PPls coming soon.
-    target : a Preliz distribution
-        This represent the prior predictive distribution **previously** elicited from the user,
+    target : a Preliz distribution or list
+        Instance of a PreliZ distribution or a list of tuples where each tuple contains a PreliZ
+        distribution and a weight.
+        This represents the prior predictive distribution **previously** elicited from the user,
         possibly using other Preliz's methods to obtain this distribution, such as maxent,
         roulette, quartile, etc.
         This should represent the domain-knowledge of the user and not any observed dataset.
-        Currently only works with a Preliz distributions. In the future we should support mixture of
-        distributions (mixture of "experts"), and maybe other options.
+    seed : int
+        A seed to initialize the Random Generator. Default is 0.
 
     Returns
     -------
@@ -50,13 +53,14 @@ def ppe(model, target):
     _log.info(""""This is an experimental method under development, use with caution.""")
 
     # Get information from PyMC model
+    rng = np.random.default_rng(seed)
     bounds, prior, p_model, var_info, var_info2, draws, free_rvs = get_model_information(model)
     # Initial point for optimization
     guess = get_guess(model, free_rvs)
     # compile PyMC model
     fmodel = compile_logp(model)
     # find prior that induce a prior predictive distribution close to target
-    prior = optimize_pymc_model(fmodel, target, draws, prior, guess, bounds, var_info, p_model)
+    prior = optimize_pymc_model(fmodel, target, draws, prior, guess, bounds, var_info, p_model, rng)
     # Fit the prior into the model's prior
     # So we can write it as a PyMC model
     new_priors = backfitting(prior, p_model, var_info2)

preliz/tests/test_ppe.py

@@ -17,19 +17,43 @@
             "sigma_z": 10,
             "target": pz.Normal(mu=174, sigma=20),
             "X": np.random.normal(0, 10, 120),
-            "new_mu_x": 174.015898,
-            "new_sigma_x": 1.971567,
-            "new_sigma_z": 19.856088,
+            "new_mu_x": 173.998287,
+            "new_sigma_x": 1.946641,
+            "new_sigma_z": 19.926345,
         },
         {
             "mu_x": [0, 1],
             "sigma_x": [10, 10],
             "sigma_z": 10,
             "target": pz.Normal(mu=174, sigma=20),
             "X": np.random.normal(0, 10, 120),
-            "new_mu_x": (174.211816, 173.822841),
-            "new_sigma_x": (2.621553, 2.691291),
-            "new_sigma_z": 19.769245,
+            "new_mu_x": (174.089111, 173.908702),
+            "new_sigma_x": (2.679288, 2.760424),
+            "new_sigma_z": 19.83069,
+        },
+        {
+            "mu_x": 0,
+            "sigma_x": 10,
+            "sigma_z": 10,
+            "target": [(pz.Normal(mu=174, sigma=20), 0.5), (pz.Normal(mu=176, sigma=19.5), 0.5)],
+            "X": np.random.normal(0, 10, 120),
+            "new_mu_x": 174.852283,
+            "new_sigma_x": 1.794118,
+            "new_sigma_z": 19.683396,
+        },
+        {
+            "mu_x": [0, 1],
+            "sigma_x": [10, 10],
+            "sigma_z": 10,
+            "target": [
+                (pz.Normal(mu=174, sigma=20), 0.5),
+                (pz.Normal(mu=176, sigma=19.5), 0.4),
+                (pz.StudentT(mu=174, sigma=20, nu=3), 0.1),
+            ],
+            "X": np.random.normal(0, 10, 120),
+            "new_mu_x": (174.75759, 174.775254),
+            "new_sigma_x": (2.728134, 2.978235),
+            "new_sigma_z": 21.247561,
         },
     ],
 )
@@ -60,10 +84,10 @@ def test_ppe(params):
             "sigma_b": 10,
             "sigma_z": 10,
             "target": pz.Normal(mu=40, sigma=7),
-            "new_mu_a": 40.013092,
-            "new_sigma_a": 0.16864,
-            "new_sigma_b": 6.991618,
-            "new_sigma_z": 6.991618,
+            "new_mu_a": 40.00908,
+            "new_sigma_a": 0.167701,
+            "new_sigma_b": 7.003267,
+            "new_sigma_z": 7.003267,
         },
     ],
 )