SciPyWrapper: Dependent and Custom Distributions¶
This notebook demonstrates independent and dependent distribution support in SciPyWrapper, including custom marginals, joint transforms, and diagnostic checks for user-defined distributions.
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
from pathlib import Path
from qmcpy import (
DigitalNetB2,
SciPyWrapper,
StudentT,
Triangular,
ZeroInflatedExpUniform,
)
from qmcpy.true_measure.triangular import TriangularDistribution
from qmcpy.util import DimensionError, ParameterError
# I like having a consistent look across all plots.
plt.rcParams["figure.figsize"] = (6, 4)
plt.rcParams["axes.grid"] = True
# Folder where all figures for the blog will be saved.
FIG_DIR = Path("figures")
FIG_DIR.mkdir(parents=True, exist_ok=True)
# ======================================================================
# 1. Acceptance rejection helpers for MC vs QMC
# ======================================================================
def sample_triangle_ar_mc(n_target, batch_size=2048, rng=None):
"""
Acceptance rejection using iid uniform proposals.
Target is the uniform distribution on the triangle:
T = {(x, y) in (0,1)^2 : 0 < y <= x < 1}.
"""
if rng is None:
rng = np.random.default_rng()
accepted = []
while sum(a.shape[0] for a in accepted) < n_target:
u = rng.random((batch_size, 2))
mask = u[:, 1] <= u[:, 0]
accepted.append(u[mask])
samples = np.concatenate(accepted, axis=0)[:n_target]
return samples
def sample_triangle_ar_qmc(n_target, batch_size=1024, seed_start=7):
"""
Acceptance rejection using QMC proposals from SciPyWrapper.
Here we wrap uniform(0,1) inside SciPyWrapper so this path uses
the same interface as all other examples.
"""
accepted = []
seed = seed_start
while sum(a.shape[0] for a in accepted) < n_target:
sampler = DigitalNetB2(2, seed=seed)
tm = SciPyWrapper(sampler, scipy_distribs=stats.uniform())
u = tm(batch_size)
mask = u[:, 1] <= u[:, 0]
accepted.append(u[mask])
seed += 1
samples = np.concatenate(accepted, axis=0)[:n_target]
return samples
# ======================================================================
# 2. Example 1: Dependent vs independent normals
# ======================================================================
def example_dependent_vs_independent_normals():
print("\n=== Example 1: dependent vs independent normals ===")
sampler = DigitalNetB2(2, seed=13)
# Independent marginals. This is essentially the original SciPyWrapper.
indep_marginals = [stats.norm(0.0, 1.0), stats.norm(0.0, 1.0)]
tm_indep = SciPyWrapper(sampler, indep_marginals)
# Joint multivariate normal with a target correlation.
rho_target = 0.8
cov = [[1.0, rho_target], [rho_target, 1.0]]
mvn = stats.multivariate_normal(mean=[0.0, 0.0], cov=cov)
tm_dep = SciPyWrapper(sampler, mvn)
n = 4096
x_indep = tm_indep(n)
x_dep = tm_dep(n)
rho_indep = np.corrcoef(x_indep.T)[0, 1]
rho_dep = np.corrcoef(x_dep.T)[0, 1]
# Functional that is very sensitive to dependence.
est_indep = np.mean(x_indep[:, 0] * x_indep[:, 1])
est_dep = np.mean(x_dep[:, 0] * x_dep[:, 1])
true_indep = 0.0 # E[X1 X2] for independent zero mean normals.
true_dep = rho_target # For MVN, E[X1 X2] = Cov(X1, X2) = rho.
print(f"Sample corr (indep) : {rho_indep:.3f}")
print(f"Sample corr (joint) : {rho_dep:.3f}")
print(f"E[X1 X2] indep est={est_indep:.3f}, true={true_indep:.3f}")
print(f"E[X1 X2] joint est={est_dep:.3f}, true={true_dep:.3f}")
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].scatter(x_indep[:, 0], x_indep[:, 1], s=6, alpha=0.4)
axes[0].set_title(f"Independent normals, rho≈{rho_indep:.2f}")
axes[0].set_xlabel("X1")
axes[0].set_ylabel("X2")
axes[1].scatter(x_dep[:, 0], x_dep[:, 1], s=6, alpha=0.4)
axes[1].set_title(f"Dependent MVN, rho≈{rho_dep:.2f}")
axes[1].set_xlabel("X1")
axes[1].set_ylabel("X2")
plt.tight_layout()
fig_path = FIG_DIR / "fig01_mvn_indep_vs_dep.png"
fig.savefig(fig_path, dpi=300, bbox_inches="tight")
plt.show()
print(f"Saved figure: {fig_path}")
# ======================================================================
# 3. Example 2: Zero inflated exponential + uniform
# ======================================================================
def example_zero_inflated_exponential():
print("\n=== Example 2: zero inflated exponential ===")
p_zero = 0.4
sampler_zi = DigitalNetB2(1, seed=21)
tm_zi = ZeroInflatedExpUniform(sampler=sampler_zi, p_zero=p_zero, lam=1.5)
n = 4096
x = tm_zi(n).ravel()
zero_rate = np.mean(x == 0.0)
print(f"Target P(X=0) : {p_zero:.3f}")
print(f"Empirical P(X=0) : {zero_rate:.3f}")
print(f"Sample mean : {x.mean():.3f}")
fig, ax = plt.subplots()
ax.hist(x, bins=60)
ax.set_title("Zero inflated exponential")
ax.set_xlabel("X")
ax.set_ylabel("count")
plt.tight_layout()
fig_path = FIG_DIR / "fig02_zero_inflated_exponential.png"
fig.savefig(fig_path, dpi=300, bbox_inches="tight")
plt.show()
print(f"Saved figure: {fig_path}")
# ======================================================================
# 4. Example 3: Acceptance rejection, MC vs QMC
# ======================================================================
def example_accept_reject_mc_vs_qmc():
print("\n=== Example 3: acceptance rejection, MC vs QMC ===")
n_target = 4096
tri_mc = sample_triangle_ar_mc(n_target)
tri_qmc = sample_triangle_ar_qmc(n_target)
rho_mc = np.corrcoef(tri_mc.T)[0, 1]
rho_qmc = np.corrcoef(tri_qmc.T)[0, 1]
print(f"corr(X, Y) MC : {rho_mc:.3f}")
print(f"corr(X, Y) QMC : {rho_qmc:.3f}")
# For uniform on the triangle T, the exact E[X Y] is 1/4.
def f(z):
return z[:, 0] * z[:, 1]
est_mc = f(tri_mc).mean()
est_qmc = f(tri_qmc).mean()
true_E = 0.25
print(f"E[X Y] MC : est={est_mc:.4f}, true={true_E:.4f}")
print(f"E[X Y] QMC : est={est_qmc:.4f}, true={true_E:.4f}")
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].scatter(tri_mc[:, 0], tri_mc[:, 1], s=6, alpha=0.4)
axes[0].set_title("Acceptance rejection with iid MC")
axes[0].set_xlabel("X")
axes[0].set_ylabel("Y")
axes[1].scatter(tri_qmc[:, 0], tri_qmc[:, 1], s=6, alpha=0.4)
axes[1].set_title("Acceptance rejection with QMC")
axes[1].set_xlabel("X")
axes[1].set_ylabel("Y")
plt.tight_layout()
fig_path = FIG_DIR / "fig03_accept_reject_mc_vs_qmc.png"
fig.savefig(fig_path, dpi=300, bbox_inches="tight")
plt.show()
print(f"Saved figure: {fig_path}")
# ======================================================================
# 5. Example 4: Custom triangular user marginal
# ======================================================================
def example_custom_triangular_marginal():
print("\n=== Example 4: custom triangular user marginal ===")
tri_custom = TriangularDistribution(c=0.3, loc=-1.0, scale=2.0) # only used for pdf overlay
sampler_tri = DigitalNetB2(1, seed=31)
tm_tri = Triangular(sampler=sampler_tri, c=0.3, loc=-1.0, scale=2.0)
n = 4096
x_tri = tm_tri(n).ravel()
print(f"Sample mean of custom triangular: {x_tri.mean():.3f}")
print(f"Sample min, max : {x_tri.min():.3f}, {x_tri.max():.3f}")
fig, ax = plt.subplots()
ax.hist(x_tri, bins=50, density=True, alpha=0.6, label="QMC samples")
xs = np.linspace(x_tri.min(), x_tri.max(), 200)
ax.plot(xs, tri_custom.pdf(xs), label="User pdf", linewidth=2.0)
ax.set_title("Custom triangular user distribution")
ax.set_xlabel("x")
ax.set_ylabel("density")
ax.legend()
plt.tight_layout()
fig_path = FIG_DIR / "fig04_custom_triangular_marginal.png"
fig.savefig(fig_path, dpi=300, bbox_inches="tight")
plt.show()
print(f"Saved figure: {fig_path}")
# ======================================================================
# 6. Example 5: Intentionally broken distribution to trigger warnings
# ======================================================================
def example_bad_custom_distribution():
print("\n=== Example 5: intentionally broken custom distribution ===")
class BadTriangular:
"""
This is a purposely broken "distribution".
- ppf is not monotone.
- pdf is constant and does not integrate to 1 on its support.
The goal is to see SciPyWrapper's sanity checks complain.
"""
def ppf(self, u):
u = np.asarray(u, dtype=float)
# Non monotone on purpose.
return np.sin(np.pi * u)
def pdf(self, x):
x = np.asarray(x, dtype=float)
# Constant density that will not normalise correctly.
return np.ones_like(x)
bad = BadTriangular()
sampler = DigitalNetB2(1, seed=99)
print("Constructing SciPyWrapper with a bad custom distribution...")
tm_bad = SciPyWrapper(sampler, bad) # warnings are expected here
# We can still draw samples, but they are not meaningful.
x = tm_bad(4)
print("Example samples from bad distribution:", x.ravel())
# ======================================================================
# 7. Example 6: Multivariate Student t (joint, dependent)
# ======================================================================
def example_multivariate_student_t_joint():
print("\n=== Example 6: multivariate Student t (joint, dependent) ===")
# Parameters
df = 5.0
rho_target = 0.8
loc = np.array([0.0, 0.0])
shape = np.array([
[1.0, rho_target],
[rho_target, 1.0],
])
# Build joint distribution + SciPyWrapper
sampler = DigitalNetB2(2, seed=41)
tm = StudentT(sampler=sampler, loc=loc, shape=shape, df=df)
n = 4096
x = tm(n)
x1 = x[:, 0]
x2 = x[:, 1]
emp_corr = np.corrcoef(x.T)[0, 1]
emp_ex12 = np.mean(x1 * x2)
# Theory (exists only if df > 2): Cov = df/(df-2) * shape
# Correlation is the same as in shape because the df/(df-2) factor cancels.
if df > 2:
true_cov12 = df / (df - 2.0) * shape[0, 1]
else:
true_cov12 = np.nan
true_corr = rho_target
print(f"df : {df:.1f}")
print(f"Target corr (from shape) : {true_corr:.3f}")
print(f"Empirical corr(X1, X2) : {emp_corr:.3f}")
print(f"E[X1 X2] empirical : {emp_ex12:.3f}")
print(f"Cov(X1, X2) theoretical : {true_cov12:.3f}")
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].hist(x1, bins=60)
axes[0].set_title("Multivariate t marginal X1")
axes[0].set_xlabel("X1")
axes[0].set_ylabel("count")
axes[1].scatter(x1, x2, s=6, alpha=0.4)
axes[1].set_title(f"Joint multivariate t, corr≈{emp_corr:.2f}")
axes[1].set_xlabel("X1")
axes[1].set_ylabel("X2")
plt.tight_layout()
fig_path = FIG_DIR / "fig05_multivariate_student_t_joint.png"
fig.savefig(fig_path, dpi=300, bbox_inches="tight")
plt.show()
print(f"Saved figure: {fig_path}")
# ======================================================================
# 8. Convenience entry point so this file runs as a script too
# ======================================================================
if __name__ == "__main__":
example_dependent_vs_independent_normals()
example_zero_inflated_exponential()
example_accept_reject_mc_vs_qmc()
example_custom_triangular_marginal()
example_bad_custom_distribution()
example_multivariate_student_t_joint()
=== Example 1: dependent vs independent normals === Sample corr (indep) : -0.000 Sample corr (joint) : 0.800 E[X1 X2] indep est=-0.000, true=0.000 E[X1 X2] joint est=0.799, true=0.800
Saved figure: figures\fig01_mvn_indep_vs_dep.png === Example 2: zero inflated exponential === Target P(X=0) : 0.400 Empirical P(X=0) : 0.400 Sample mean : 0.400
C:\Users\Owner\Downloads\QMCSoftware\qmcpy\true_measure\scipy_wrapper.py:341: UserWarning: Custom univariate distribution has no 'pdf' or 'logpdf'. Weights will be treated as 1 for this marginal. warnings.warn(
Saved figure: figures\fig02_zero_inflated_exponential.png === Example 3: acceptance rejection, MC vs QMC === corr(X, Y) MC : 0.514 corr(X, Y) QMC : 0.501 E[X Y] MC : est=0.2554, true=0.2500 E[X Y] QMC : est=0.2505, true=0.2500
Saved figure: figures\fig03_accept_reject_mc_vs_qmc.png === Example 4: custom triangular user marginal === Sample mean of custom triangular: -0.133 Sample min, max : -0.993, 0.976
Saved figure: figures\fig04_custom_triangular_marginal.png === Example 5: intentionally broken custom distribution === Constructing SciPyWrapper with a bad custom distribution... Example samples from bad distribution: [0.99779416 0.30984091 0.43378429 0.77064796] === Example 6: multivariate Student t (joint, dependent) === df : 5.0 Target corr (from shape) : 0.800 Empirical corr(X1, X2) : 0.808 E[X1 X2] empirical : 1.389 Cov(X1, X2) theoretical : 1.333
C:\Users\Owner\Downloads\QMCSoftware\qmcpy\true_measure\scipy_wrapper.py:228: UserWarning: SciPyWrapper received a custom univariate distribution (not a scipy.stats frozen distribution) and it failed sanity checks: - ppf() is not nondecreasing (looks non-monotone) self._setup_marginals(scipy_distribs)
Saved figure: figures\fig05_multivariate_student_t_joint.png
Interpreting the above outputs¶
Example 1: Dependent vs Independent Normals¶
This example checks whether SciPyWrapper correctly handles both independent marginal distributions and a dependent multivariate normal distribution.
In the independent case, the two variables are sampled separately:
$$ X_1 \sim N(0,1), \qquad X_2 \sim N(0,1). $$
Since they are independent, we expect
$$ \operatorname{Corr}(X_1, X_2) \approx 0 $$
and
$$ \mathbb{E}[X_1X_2] = \mathbb{E}[X_1]\mathbb{E}[X_2] = 0. $$
The sample correlation and the estimate of $\mathbb{E}[X_1X_2]$ being close to zero confirms that the independent marginal case is behaving as expected.
In the dependent case, the samples come from a multivariate normal distribution with covariance matrix
$$ \Sigma = \begin{pmatrix} 1 & \rho \ \rho & 1 \end{pmatrix}. $$
Here, $\rho$ controls the dependence between $X_1$ and $X_2$. Since both variables have mean zero,
$$ \mathbb{E}[X_1X_2] = \operatorname{Cov}(X_1, X_2) = \rho. $$
The dependent sample correlation and the estimate of $\mathbb{E}[X_1X_2]$ are close to the target value $\rho$, which shows that the wrapper is preserving the intended dependence. The diagonal pattern in the scatter plot also visually confirms the positive correlation.
Example 2: Zero-Inflated Exponential¶
This example checks the one-dimensional zero-inflated exponential transform. The distribution has a point mass at zero and an exponential tail for positive values.
Let
$$ U \sim \operatorname{Uniform}(0,1), $$
with zero-inflation probability $p_0$ and exponential rate $\lambda$. The transform is
$$ X = \begin{cases} 0, & U \le p_0, \ -\frac{1}{\lambda}\log\left(1-\frac{U-p_0}{1-p_0}\right), & U > p_0. \end{cases} $$
This gives
$$ \mathbb{P}(X=0)=p_0. $$
For this example, $p_0=0.4$, so about 40% of the samples should be exactly zero. The large spike at zero in the histogram confirms the point mass.
The remaining samples come from the exponential part of the distribution, so the histogram should have a decreasing right tail. The theoretical mean is
$$ \mathbb{E}[X] = (1-p_0)\frac{1}{\lambda}. $$
With $p_0=0.4$ and $\lambda=1.5$,
$$ \mathbb{E}[X] = \frac{0.6}{1.5} = 0.4. $$
The empirical zero rate and sample mean both match the theoretical values closely, which shows that the one-dimensional transform is producing the intended zero-inflated exponential distribution.
Example 3: Acceptance-Rejection Sampling, MC vs QMC¶
This example compares ordinary Monte Carlo acceptance-rejection sampling with a QMC-based acceptance-rejection approach.
The target distribution is uniform on the triangular region
$$ T = {(x,y): 0 < y \le x < 1}. $$
The acceptance rule keeps only the points satisfying
$$ y \le x. $$
Therefore, both the MC and QMC scatter plots should fill the same triangular region. This confirms that the acceptance-rejection rule is selecting points from the correct target domain.
The example also estimates
$$ \mathbb{E}[XY] $$
for a uniform distribution on this triangle. The exact value is
$$ \mathbb{E}[XY] = \frac{1}{4}. $$
Both estimates being close to $0.25$ confirms that the accepted samples are consistent with the target triangular distribution. The QMC estimate may often look more stable because QMC points are more evenly distributed before the rejection step.
The positive correlation in the accepted samples is also expected. Since accepted points must satisfy $y \le x$, larger values of $x$ allow larger possible values of $y$, creating positive dependence between $X$ and $Y$.
Example 4: Custom Triangular Marginal¶
This example checks whether SciPyWrapper can use a custom one-dimensional distribution that provides a ppf and pdf.
The triangular distribution is supported on
$$ [a,b] = [\text{loc}, \text{loc}+\text{scale}], $$
with mode
$$ m = \text{loc} + c \cdot \text{scale}. $$
The inverse CDF maps uniform samples from $[0,1]$ into samples from the triangular distribution. The histogram should match the analytic PDF: it increases up to the mode and decreases after the mode.
The mean of a triangular distribution is
$$ \mathbb{E}[X] = \frac{a+b+m}{3}. $$
In this example,
$$ a=-1, \qquad b=1, \qquad m=-1+0.3(2)=-0.4. $$
Therefore,
$$ \mathbb{E}[X] = \frac{-1+1-0.4}{3} \approx -0.133. $$
The samples stay inside the expected interval, and the sample mean is close to the theoretical mean. The histogram also follows the analytic PDF, showing that the custom ppf and pdf are consistent.
Example 5: Intentionally Broken Custom Distribution¶
This example intentionally passes a bad custom distribution into SciPyWrapper to test whether the sanity checks catch invalid behavior.
The broken distribution uses a non-monotone inverse CDF:
$$ \operatorname{ppf}(u)=\sin(\pi u). $$
A valid inverse CDF should be nondecreasing, because larger probability values should not map backward to smaller sample values. This distribution also uses a constant PDF that is not properly normalized for the implied distribution.
Warnings are expected in this example. The purpose is not to generate meaningful samples. Instead, this example demonstrates that SciPyWrapper can detect suspicious custom distributions and warn the user.
The notebook continuing to run after showing warnings is the desired behavior here. The warning messages indicate that the sanity checks are catching the intentionally invalid distribution.
Example 6: Multivariate Student-t Distribution¶
This example checks whether the dependent multivariate Student-t transform is working correctly.
The multivariate Student-t distribution uses a shape matrix to represent dependence:
$$ \Sigma = \begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix}. $$
The parameter $\rho$ controls the target correlation between the two coordinates. Therefore, the empirical correlation from the generated samples should be close to $\rho$.
For a Student-t distribution with degrees of freedom $\nu > 2$, the covariance matrix is
$$ \operatorname{Cov}(X) = \frac{\nu}{\nu-2}\Sigma. $$
So the off-diagonal covariance is
$$ \operatorname{Cov}(X_1,X_2)=\frac{\nu}{\nu-2}\rho. $$
With $\nu=5$ and $\rho=0.8$,
$$ \operatorname{Cov}(X_1,X_2)=\frac{5}{3}(0.8)\approx 1.333. $$
Since the variables are centered at zero,
$$ \mathbb{E}[X_1X_2] = \operatorname{Cov}(X_1,X_2). $$
The empirical correlation being close to the target $\rho$ confirms that the Student-t transform is preserving the intended dependence. The scatter plot shows this dependence visually, while the marginal histogram shows the heavier-tailed behavior expected from a Student-t distribution. Some variation in the covariance estimate is natural because Student-t samples can include more extreme values.
def example1_zero_inflated_exponential():
print("\n=== Example 2: zero-inflated exponential ===")
p_zero = 0.4
sampler_zi = DigitalNetB2(1, replications=3, seed=21)
tm_zi = ZeroInflatedExpUniform(
sampler=sampler_zi,
p_zero=p_zero,
lam=1.5,
)
n = 8
u = sampler_zi(n)
np.set_printoptions(precision=6, suppress=True)
#print(u)
print("Shape of u:", u.shape)
print(u.squeeze(-1))
# Shape: (replications, samples per replication, dimension) = (r, n, 1)
x = tm_zi(n)
print("Shape of x:", x.shape)
zero_rate_by_replication = np.mean(x == 0.0, axis=(1, 2))
mean_by_replication = np.mean(x, axis=(1, 2))
# Combine all replications only for one overall histogram and summary.
x_flat = x.reshape(-1)
print(f"Target P(X=0) : {p_zero:.3f}")
print(f"Overall empirical P(X=0) : {np.mean(x_flat == 0.0):.3f}")
print(f"Overall sample mean : {x_flat.mean():.3f}")
fig, ax = plt.subplots()
ax.hist(x_flat, bins=60)
ax.set_title("Zero-inflated exponential")
ax.set_xlabel("X")
ax.set_ylabel("count")
plt.tight_layout()
plt.show()
if __name__ == "__main__":
example1_zero_inflated_exponential()
=== Example 2: zero-inflated exponential === Shape of u: (3, 8, 1) [[0.494023 0.826311 0.076979 0.74325 0.29804 0.881792 0.130873 0.548609] [0.562424 0.216008 0.795917 0.451026 0.739626 0.023183 0.974154 0.257162] [0.462443 0.708298 0.092731 0.83027 0.314138 0.575607 0.194165 0.947348]] Shape of x: (3, 8, 1) Target P(X=0) : 0.400 Overall empirical P(X=0) : 0.375 Overall sample mean : 0.403
