Add support for mixed distributions in TrueMeasure

[Laasya-73]

The goal is to support a mixed distribution composed of several component distributions, each selected with a specified probability.

Suppose there are $s$ component distributions with densities

$$\rho_0(\boldsymbol{x}), \rho_1(\boldsymbol{x}), \ldots, \rho_{s-1}(\boldsymbol{x}), \qquad \boldsymbol{x} \in \mathbb{R}^d,$$

and corresponding probabilities

$$p_0, p_1, \ldots, p_{s-1},$$

satisfying

$$p_j > 0, \qquad j = 0,\ldots,s-1,$$

and

$$\sum_{j=0}^{s-1} p_j = 1.$$

The density of the resulting mixed distribution is

$$\rho(\boldsymbol{x}) = \sum_{j=0}^{s-1} p_j \rho_j(\boldsymbol{x}).$$

Equivalently,

$$\rho(\boldsymbol{x}) = p_0\rho_0(\boldsymbol{x}) + p_1\rho_1(\boldsymbol{x}) + \cdots + p_{s-1}\rho_{s-1}(\boldsymbol{x}).$$

Proposed sampling transform

Sampling from the mixture requires one additional uniform coordinate.

Let

$$\boldsymbol{u} = (u_0,u_1,\ldots,u_d) \in [0,1]^{d+1}.$$

The first coordinate, $u_0$, is used to determine which component distribution is selected.

Define the cumulative probabilities

$$P_j = \sum_{k=0}^{j} p_k, \qquad j=0,\ldots,s-1,$$

with

$$P_{-1}=0$$

and

$$P_{s-1}=1.$$

Component $j$ is selected when

$$P_{j-1} \le u_0 < P_j,$$

for

$$j=0,\ldots,s-2.$$

For the final component, $j=s-1$, use

$$P_{s-2} \le u_0 \le P_{s-1}=1.$$

This ensures that the full interval

$$u_0 \in [0,1]$$

is covered, including the boundary values $u_0=0$ and $u_0=1$.

Once component $j$ has been selected, define

$$\boldsymbol{v} = (u_1,\ldots,u_d) \in [0,1]^d.$$

The remaining $d$ coordinates are transformed using the selected component distribution.

If the transform associated with component $j$ is

$$T_j: [0,1]^d \rightarrow \mathbb{R}^d,$$

then the final sample is

$$\boldsymbol{x} = T_j(\boldsymbol{v}) = T_j(u_1,\ldots,u_d).$$

Therefore, the complete mixed-distribution transform is

$$T: [0,1]^{d+1} \rightarrow \mathbb{R}^d.$$

The first coordinate determines the selected component,

$$u_0 \longrightarrow j,$$

while the remaining coordinates generate the sample from that component,

$$(u_1,\ldots,u_d) \longrightarrow T_j(u_1,\ldots,u_d).$$

Initial implementation

The initial implementation should support:

• Multiple component distributions
• User-provided probabilities $p_0,\ldots,p_{s-1}$
• Validation that all probabilities are positive:

$$p_j > 0, \qquad j=0,\ldots,s-1$$

• Validation that the probabilities sum to $1$:

$$\sum_{j=0}^{s-1} p_j = 1$$

• Validation that all component distributions have compatible output dimension $d$
• Component selection using the additional coordinate $u_0$
• Transformation of $(u_1,\ldots,u_d)$ through the selected component
• Evaluation of the mixed density
• Corresponding TrueMeasure weight calculation
• Tests for expected sampling behavior
• Tests for invalid probabilities
• Tests for incompatible component dimensions
• Tests for boundary cases in component selection

Dimension handling

An important implementation detail is that the input dimension and output dimension are different.

The underlying sampler generates points in

$$[0,1]^{d+1},$$

so the sampler dimension is

$$d+1.$$

However, the resulting mixed-distribution samples belong to

$$\mathbb{R}^d,$$

so the TrueMeasure output dimension is

$$d.$$

Conceptually,

(u_0, u_1, ..., u_d)
  |      |
  |      +---- (u_1, ..., u_d) --> component transform T_j
  |
  +----------- u_0 --> select component j

                           |
                           v

                  x = T_j(u_1, ..., u_d)

                           |
                           v

                         R^d

Thus,

$$\text{sampler dimension}=d+1,$$

while

$$\text{TrueMeasure output dimension}=d.$$

This may require special handling because AbstractTrueMeasure currently derives its dimension directly from the sampler.

The mixed TrueMeasure therefore needs to distinguish between:

• the sampler dimension, $d+1$, and
• the transformed sample dimension, $d$.

Mixed density

For any

$$\boldsymbol{x}\in\mathbb{R}^d,$$

the density of the mixture is

$$\rho(\boldsymbol{x}) = \sum_{j=0}^{s-1} p_j\rho_j(\boldsymbol{x}).$$

The density must include contributions from all component distributions, not only the component that was selected when generating a particular sample.

Therefore, the corresponding TrueMeasure weight calculation should use the complete mixture density

$$\rho(\boldsymbol{x}),$$

rather than only the density of the selected component,

$$\rho_j(\boldsymbol{x}).$$

[nomadbl]

The initial implementation outlined above is not a proper QMC procedure, but rather a mixture of MC and QMC. The results would be wrong because removing samples from the individual components (due to the sub-sampling outlined above) would break the symmetry of the lattice or other grid used. i.e. it would be a sum of incorrect integrals.

This proposed feature can be alternatively approached by:

1. computing the integrals corresponding to the individual components
2. summing with the component weights

The complication might be how to implement this in the architecture of qmcpy.

[fjhickernell]

The initial implementation outlined above is not a proper QMC procedure, but rather a mixture of MC and QMC. The results would be wrong because removing samples from the individual components (due to the sub-sampling outlined above) would break the symmetry of the lattice or other grid used. i.e. it would be a sum of incorrect integrals.

This proposed feature can be alternatively approached by:

1. computing the integrals corresponding to the individual components
2. summing with the component weights

The complication might be how to implement this in the architecture of qmcpy.

Thanks for your comment, Lior. You have a valid point. We will look into this and we reply.