> ## Documentation Index
> Fetch the complete documentation index at: https://docs.haiqu.ai/llms.txt
> Use this file to discover all available pages before exploring further.

# Multivariate Probability Distributions

> Load multidimensional probability densities and copula-based joint probability density functions into quantum states.

## Multivariate Distribution Loading

Haiqu SDK supports loading named multidimensional probability density functions (PDFs) into quantum circuits with `haiqu.multivariate_distribution_loading`.

```python theme={null}
job = haiqu.multivariate_distribution_loading(
    distribution_name="bivariate_student_t",
    num_qubits=8,  # 4 qubits for x and 4 qubits for y
    interval=((-4.0, 4.0), (-4.0, 4.0)),
    encoding="probability",
    distribution_params={
        "mu": [0.0, 0.0],
        "sigma": [[1.0, 0.3], [0.3, 1.0]],
        "nu": 5.0,
    },
)

data_loading_gate = job.result()  # a HaiquCircuitGate, encapsulating state preparation circuit
```

<Tip>
  Circuit fidelity is stored in `job.quality`. This is the MPS-to-circuit compilation fidelity.
</Tip>

For a $d$-dimensional distribution with $n_j$ qubits allocated to dimension $j$, the qubit layout groups dimensions sequentially: dimension $1$ occupies qubits $q^{(1)}_0, \ldots, q^{(1)}_{n_1-1}$, dimension $2$ occupies $q^{(2)}_0, \ldots, q^{(2)}_{n_2-1}$, and so on up to dimension $d$, where $q^{(j)}_0$ is the least-significant-bit qubit for dimension $j$. If `num_qubits=(n_1, ..., n_d)`, the circuit uses $n_j$ qubits for dimension $j$.

Each qubit group encodes a grid index in binary: for dimension $j$, the integer index

$i_j = q^{(j)}_0 \cdot 2^0 + q^{(j)}_1 \cdot 2^1 + \cdots + q^{(j)}_{n_j-1} \cdot 2^{n_j-1}$

maps to the grid point

$x^{(j)}_{i_j} = \frac{\mathrm{high}_j - \mathrm{low}_j}{2^{n_j} - 1} i_j + \mathrm{low}_j \in [\mathrm{low}_j, \mathrm{high}_j].$

Qiskit uses the qubit-ordering convention $|q_{N-1}, \ldots, q_1, q_0\rangle$ (most-significant qubit on the left), so the statevector index ordering is the reverse of the layout above: the full register reads $|q^{(d)}_{n_d-1}, \ldots, q^{(d)}_0, \ldots, q^{(1)}_{n_1-1}, \ldots, q^{(1)}_0\rangle$ in Qiskit notation.

## Encoding Modes

| Mode                     | Amplitudes                                           | Measurement probabilities                 | Use case                          |
| :----------------------- | :--------------------------------------------------- | :---------------------------------------- | :-------------------------------- |
| `encoding="amplitude"`   | proportional to the density value                    | proportional to the squared density value | Direct function-amplitude loading |
| `encoding="probability"` | proportional to the square root of the density value | proportional to the density value         | Probability distribution loading  |

For a $d$-dimensional grid point $\mathbf{x}_\mathbf{i} = (x^{(1)}_{i_1}, \ldots, x^{(d)}_{i_d})$ and density $f$, the two modes encode:

$\psi_{\mathbf{i}}^{\mathrm{amplitude}} = \frac{f(\mathbf{x}_\mathbf{i})}{\sqrt{\sum_{\mathbf{k}} f(\mathbf{x}_\mathbf{k})^2}}$

$\psi_{\mathbf{i}}^{\mathrm{probability}} = \sqrt{\frac{f(\mathbf{x}_\mathbf{i})}{\sum_{\mathbf{k}} f(\mathbf{x}_\mathbf{k})}}$

<Tip>
  Use `encoding="probability"` when `distribution_name` refers to a PDF and measurement outcomes should sample from that PDF.
</Tip>

## Parameters

| Parameter                      | Type                                        | Description                                                                                   |
| :----------------------------- | :------------------------------------------ | :-------------------------------------------------------------------------------------------- |
| `distribution_name`            | `str`                                       | Name of a supported distribution or copula.                                                   |
| `num_qubits`                   | `int` or tuple of `int`                     | Total qubit count split evenly across dimensions, or a per-dimension tuple `(n_1, ..., n_d)`. |
| `interval`                     | `(a, b)` or `((a_1, b_1), ..., (a_d, b_d))` | Shared interval for all dimensions, or per-dimension intervals. Defaults to `(0, 1)`.         |
| `encoding`                     | `str`                                       | `"amplitude"` or `"probability"`. Defaults to `"amplitude"`.                                  |
| `distribution_params`          | `dict`                                      | Parameters for direct 2D distributions.                                                       |
| `copula_params`                | `dict`                                      | Parameters for a named copula.                                                                |
| `marginal_distribution_names`  | sequence of `str`                           | Two `scipy.stats` marginal names for copula-based joint PDFs.                                 |
| `marginal_distribution_params` | sequence of `dict`                          | Parameters for each marginal distribution.                                                    |
| `num_layers`                   | `int`                                       | Number of disentangler layers in MPS-to-circuit compilation. Defaults to `1`.                 |
| `truncation_cutoff`            | `float` or `None`                           | Entanglement cutoff for later layers. Defaults to `1e-6`.                                     |
| `fine_tuning_iterations`       | `int`                                       | Maximum gradient-based refinement iterations. Defaults to `20`.                               |

## Supported Direct 2D Distributions

| Name                              | Support                                          | Parameters                                                 |
| :-------------------------------- | :----------------------------------------------- | :--------------------------------------------------------- |
| `bivariate_student_t`             | all real `x` and `y`                             | `mu`, `sigma`, `nu`                                        |
| `bivariate_gamma`                 | positive `x` and `y`                             | `nu`, `rho`                                                |
| `bivariate_von_mises`             | angular domain, usually `[-pi, pi) x [-pi, pi)`  | `mu1`, `mu2`, `kappa1`, `kappa2`, `lam`, optional `n_norm` |
| `marshall_olkin_weibull`          | positive `x` and `y`                             | `lambda1`, `lambda2`, `lambda3`, `k`                       |
| `bivariate_poisson_normal_approx` | continuous approximation on all real `x` and `y` | `lambda1`, `lambda2`, `lambda3`                            |

### Bivariate Student-t

**Name in SDK:** `bivariate_student_t`

A heavy-tail generalization of the bivariate Gaussian. It is useful when the joint distribution needs Gaussian-like correlation structure while retaining robustness to tail events and outliers.

The density is the standard two-dimensional Student-t PDF with location vector `mu`, scale matrix `sigma`, and degrees of freedom `nu`.

$f(x,y) = \frac{\Gamma((\nu+2)/2)}{\Gamma(\nu/2)\nu\pi|\Sigma|^{1/2}}\left(1 + \frac{(z-\mu)^T\Sigma^{-1}(z-\mu)}{\nu}\right)^{-(\nu+2)/2}$

where $z=(x,y)^T$.

| Parameter | Type                      | Constraints       | Description        |
| :-------- | :------------------------ | :---------------- | :----------------- |
| `mu`      | array-like shape `(2,)`   | finite            | Location vector    |
| `sigma`   | array-like shape `(2, 2)` | positive-definite | Scale matrix       |
| `nu`      | `float`                   | greater than `0`  | Degrees of freedom |

```python theme={null}
data_loading_gate = haiqu.multivariate_distribution_loading(
    distribution_name="bivariate_student_t",
    num_qubits=(4, 4),
    interval=((-4.0, 4.0), (-4.0, 4.0)),
    encoding="probability",
    distribution_params={
        "mu": [0.0, 0.0],
        "sigma": [[1.0, 0.5], [0.5, 1.0]],
        "nu": 4.0,
    },
).result()
```

### Bivariate Gamma

**Name in SDK:** `bivariate_gamma`

Kibble's bivariate Gamma distribution with Gamma marginals and a positive dependence parameter. The implementation uses the scaled modified Bessel function for numerical stability.

$f(x,y) = \frac{(xy)^{(\nu-1)/2}}{\Gamma(\nu)(1-\rho)\rho^{(\nu-1)/2}}\exp\left(-\frac{x+y}{1-\rho}\right)I_{\nu-1}\left(\frac{2\sqrt{\rho xy}}{1-\rho}\right)$

where $I_{\nu-1}$ is the modified Bessel function of the first kind. For `rho=0`, the implementation returns the product of two independent Gamma densities.

| Parameter | Type    | Constraints      | Description                    |
| :-------- | :------ | :--------------- | :----------------------------- |
| `nu`      | `float` | greater than `0` | Marginal Gamma shape           |
| `rho`     | `float` | in `[0, 1)`      | Positive correlation parameter |

```python theme={null}
data_loading_gate = haiqu.multivariate_distribution_loading(
    distribution_name="bivariate_gamma",
    num_qubits=8,
    interval=((0.0, 8.0), (0.0, 8.0)),
    encoding="probability",
    distribution_params={"nu": 2.0, "rho": 0.4},
).result()
```

### Bivariate von Mises

**Name in SDK:** `bivariate_von_mises`

A circular analogue of the bivariate Gaussian, defined on angular variables. The implementation uses the Sine model and computes the normalizing constant numerically.

$f(\theta_1,\theta_2) = \exp\left(\kappa_1\cos(\theta_1-\mu_1) + \kappa_2\cos(\theta_2-\mu_2) + \lambda\sin(\theta_1-\mu_1)\sin(\theta_2-\mu_2)\right)$

| Parameter | Type    | Constraints  | Description                                                        |
| :-------- | :------ | :----------- | :----------------------------------------------------------------- |
| `mu1`     | `float` | radians      | Mean direction of the first angle                                  |
| `mu2`     | `float` | radians      | Mean direction of the second angle                                 |
| `kappa1`  | `float` | at least `0` | Concentration of the first angle                                   |
| `kappa2`  | `float` | at least `0` | Concentration of the second angle                                  |
| `lam`     | `float` | real         | Coupling parameter                                                 |
| `n_norm`  | `int`   | positive     | Grid size per dimension for numerical normalization, default `100` |

```python theme={null}
data_loading_gate = haiqu.multivariate_distribution_loading(
    distribution_name="bivariate_von_mises",
    num_qubits=8,
    interval=((-3.14159, 3.14159), (-3.14159, 3.14159)),
    encoding="probability",
    distribution_params={
        "mu1": 0.0,
        "mu2": 0.4,
        "kappa1": 2.0,
        "kappa2": 1.5,
        "lam": 0.8,
    },
).result()
```

### Marshall-Olkin Bivariate Weibull

**Name in SDK:** `marshall_olkin_weibull`

A shock model with independent failure sources for each variable and a shared common shock. The implemented density is the absolutely continuous part; the singular mass on the diagonal `x = y` is not represented.

For `x < y`, the density is:

$f(x,y) = \lambda_1(\lambda_2+\lambda_3)k^2(xy)^{k-1}\exp\left(-\lambda_1x^k-(\lambda_2+\lambda_3)y^k\right)$

For `x > y`, the density is:

$f(x,y) = \lambda_2(\lambda_1+\lambda_3)k^2(xy)^{k-1}\exp\left(-(\lambda_1+\lambda_3)x^k-\lambda_2y^k\right)$

On `x = y`, the absolutely continuous density returned by this function is `0`.

| Parameter | Type    | Constraints      | Description                    |
| :-------- | :------ | :--------------- | :----------------------------- |
| `lambda1` | `float` | greater than `0` | Independent shock rate for `X` |
| `lambda2` | `float` | greater than `0` | Independent shock rate for `Y` |
| `lambda3` | `float` | at least `0`     | Common-shock rate              |
| `k`       | `float` | greater than `0` | Weibull shape                  |

```python theme={null}
data_loading_gate = haiqu.multivariate_distribution_loading(
    distribution_name="marshall_olkin_weibull",
    num_qubits=8,
    interval=((0.0, 4.0), (0.0, 4.0)),
    encoding="probability",
    distribution_params={
        "lambda1": 1.0,
        "lambda2": 1.2,
        "lambda3": 0.4,
        "k": 1.5,
    },
).result()
```

### Bivariate Poisson Normal Approximation

**Name in SDK:** `bivariate_poisson_normal_approx`

A continuous Gaussian approximation to a bivariate Poisson model constructed from `X = Z_1 + Z_3` and `Y = Z_2 + Z_3`, where each `Z_i` is Poisson-distributed. The shared component `Z_3` induces covariance.

$\mu = [\lambda_1+\lambda_3,\ \lambda_2+\lambda_3]$

$\Sigma_{11} = \lambda_1+\lambda_3,\quad \Sigma_{22} = \lambda_2+\lambda_3,\quad \Sigma_{12} = \Sigma_{21} = \lambda_3$

| Parameter | Type    | Constraints      | Description                         |
| :-------- | :------ | :--------------- | :---------------------------------- |
| `lambda1` | `float` | greater than `0` | Rate of `X`'s independent component |
| `lambda2` | `float` | greater than `0` | Rate of `Y`'s independent component |
| `lambda3` | `float` | at least `0`     | Shared component rate               |

```python theme={null}
data_loading_gate = haiqu.multivariate_distribution_loading(
    distribution_name="bivariate_poisson_normal_approx",
    num_qubits=8,
    interval=((0.0, 12.0), (0.0, 12.0)),
    encoding="probability",
    distribution_params={
        "lambda1": 2.0,
        "lambda2": 3.0,
        "lambda3": 1.0,
    },
).result()
```

## Copula-Based Joint PDFs

Copula loading builds a joint PDF from two `scipy.stats` marginals and a supported [copula](https://en.wikipedia.org/wiki/Copula_\(statistics\)) density using [Sklar's theorem](https://en.wikipedia.org/wiki/Copula_\(statistics\)#Sklar's_theorem).

$f(x,y) = c(F_X(x), F_Y(y)) f_X(x) f_Y(y)$

Use a copula name in `distribution_name`, provide the copula parameters in `copula_params`, and provide exactly two marginal names and parameter dictionaries.

| Copula name       | Parameters           | Dependence behavior                             |
| :---------------- | :------------------- | :---------------------------------------------- |
| `gaussian_copula` | `rho` in `(-1, 1)`   | Gaussian dependence structure                   |
| `clayton_copula`  | `theta` at least `0` | Lower-tail dependence                           |
| `gumbel_copula`   | `theta` at least `1` | Upper-tail dependence                           |
| `frank_copula`    | non-zero `theta`     | Symmetric dependence without tail concentration |

The supported copula densities are:

$c_{\mathrm{Gaussian}}(u,v) = \frac{1}{\sqrt{1-\rho^2}}\exp\left(\frac{\rho(2xy-\rho(x^2+y^2))}{2(1-\rho^2)}\right),\quad x=\Phi^{-1}(u),\ y=\Phi^{-1}(v)$

$c_{\mathrm{Clayton}}(u,v) = (1+\theta)(uv)^{-(\theta+1)}(u^{-\theta}+v^{-\theta}-1)^{-1/\theta-2}$

$c_{\mathrm{Gumbel}}(u,v) = \exp(-A^{1/\theta})(st)^{\theta-1}A^{1/\theta-2}\frac{A^{1/\theta}+\theta-1}{uv},\quad s=-\ln u,\ t=-\ln v,\ A=s^\theta+t^\theta$

$c_{\mathrm{Frank}}(u,v) = \frac{-\theta(e^{-\theta}-1)e^{-\theta(u+v)}}{\left(e^{-\theta}-1+(e^{-\theta u}-1)(e^{-\theta v}-1)\right)^2}$

```python theme={null}
data_loading_gate = haiqu.multivariate_distribution_loading(
    distribution_name="gaussian_copula",
    num_qubits=(4, 4),
    interval=((-4.0, 4.0), (0.0, 8.0)),
    encoding="probability",
    copula_params={"rho": 0.6},
    marginal_distribution_names=["norm", "gamma"],
    marginal_distribution_params=[
        {"loc": 0.0, "scale": 1.0},
        {"a": 2.0, "scale": 1.0},
    ],
).result()
```

<Note>
  `distribution_params` are only valid for direct 2D distributions. For copulas, use `copula_params`, `marginal_distribution_names`, and `marginal_distribution_params`.
</Note>

## Normal and Lognormal Distributions

Normal and lognormal marginals can be loaded via the copula interface using `scipy.stats` distribution names. Both examples below use a Gaussian copula, which preserves the familiar Gaussian dependence structure.

### Bivariate Normal

Use `"norm"` as both marginal names. The `loc` and `scale` parameters correspond to the mean and standard deviation of each margin.

```python theme={null}
job = haiqu.multivariate_distribution_loading(
    distribution_name="gaussian_copula",
    num_qubits=(4, 4),
    interval=((-4.0, 4.0), (-4.0, 4.0)),
    encoding="probability",
    copula_params={"rho": 0.0},
    marginal_distribution_names=["norm", "norm"],
    marginal_distribution_params=[
        {"loc": 0.0, "scale": 1.0},
        {"loc": 0.0, "scale": 1.0},
    ],
)

data_loading_gate = job.result()
```

### Bivariate Lognormal

Use `"lognorm"` as both marginal names. In `scipy.stats`, `lognorm` is parameterized by shape `s` (the standard deviation of the underlying normal in log-space), `loc` (shift, usually `0`), and `scale` (equal to `exp(mu)` where `mu` is the mean of the underlying normal).

```python theme={null}
job = haiqu.multivariate_distribution_loading(
    distribution_name="gaussian_copula",
    num_qubits=(4, 4),
    interval=((0.01, 6.0), (0.01, 6.0)),
    encoding="probability",
    copula_params={"rho": 0.4},
    marginal_distribution_names=["lognorm", "lognorm"],
    marginal_distribution_params=[
        {"s": 0.5, "loc": 0.0, "scale": 1.0},
        {"s": 0.8, "loc": 0.0, "scale": 1.0},
    ],
)

data_loading_gate = job.result()
```

<Note>
  The lognormal interval must exclude zero. Use a small positive lower bound such as `0.01` rather than `0.0` to avoid evaluating the density at the boundary.
</Note>