API Reference¶

This page provides comprehensive documentation for all bezierv classes, functions, and algorithms.

Core Classes¶

Bezierv - Main Distribution Class¶

The primary class for representing and working with Bézier random variables.

bezierv.classes.bezierv.Bezierv ¶

Bezierv(n: int, controls_x: ArrayLike = None, controls_z: ArrayLike = None)

Initialize a Bézier random variable.

Sets up control points, binomial coefficient arrays, and consecutive difference arrays. Moment attributes are set to np.inf as a sentinel indicating they have not yet been computed.

Parameters:

Name	Type	Description	Default
`n`	`int`	Degree of the Bézier curve. The curve has `n + 1` control points (paper notation: degree-`n` Bézier distribution).	required
`controls_x`	`array-like of shape (n + 1,)`	x-coordinates of the control points. Must be strictly increasing. If `None`, a zero array is created (placeholder state).	`None`
`controls_z`	`array-like of shape (n + 1,)`	z-coordinates (CDF values) of the control points. Must be non-decreasing with `z[0] = 0` and `z[n] = 1`. If `None`, a zero array is created (placeholder state).	`None`

Attributes:

Name	Type	Description
`n`	`int`	Degree of the Bézier curve.
`deltas_x`	`(ndarray, shape(n))`	Differences between consecutive x control points.
`deltas_z`	`(ndarray, shape(n))`	Differences between consecutive z control points.
`comb`	`(ndarray, shape(n + 1))`	Binomial coefficients `C(n, i)` for `i = 0, …, n`.
`comb_minus`	`(ndarray, shape(n))`	Binomial coefficients `C(n-1, i)` for `i = 0, …, n-1` (used for PDF evaluation).
`support`	`tuple of float`	`(controls_x[0], controls_x[-1])`, or `(-inf, inf)` in placeholder state.
`controls_x`	`(ndarray, shape(n + 1))`	x-coordinates of the control points.
`controls_z`	`(ndarray, shape(n + 1))`	z-coordinates of the control points.
`raw_moments`	`dict[int, float]`	Cache of computed raw moments `E[X^r]` keyed by order `r`. Cleared by `update_bezierv`.

Methods:¶

update_bezierv ¶

update_bezierv(controls_x: array, controls_z: array)

Update the control points, support, and delta arrays in-place.

Parameters:

Name	Type	Description	Default
`controls_x`	`array-like of shape (n + 1,)`	New x-coordinates of the control points. Must be strictly increasing.	required
`controls_z`	`array-like of shape (n + 1,)`	New z-coordinates of the control points. Must be non-decreasing.	required

combinations ¶

combinations()

Compute and store binomial coefficients.

deltas ¶

deltas()

Compute the differences between consecutive control points.

bernstein ¶

bernstein(t: float, i: int, combinations: ndarray, n: int) -> float

Evaluate the i-th degree-n Bernstein basis polynomial at t.

Computes B_{n,i}(t) = C(n,i) * t^i * (1-t)^(n-i) as defined in Eq. (1) of the paper.

Parameters:

Name	Type	Description	Default
`t`	`float`	Parameter value in `[0, 1]`.	required
`i`	`int`	Index of the basis polynomial, `0 <= i <= n`.	required
`combinations`	`ndarray`	Precomputed binomial coefficients; pass `self.comb` for degree `n` or `self.comb_minus` for degree `n - 1`.	required
`n`	`int`	Degree of the Bernstein polynomial.	required

Returns:

Type	Description
`float`	Value of `B_{n,i}(t)`.

poly_x ¶

poly_x(t: float, controls_x: ndarray = None) -> float

Evaluate the x-component of the Bézier curve at parameter t.

Computes B_x(t) = ∑ᵢ B_{n,i}(t) * xᵢ (Eq. 2 of the paper).

Parameters:

Name	Type	Description	Default
`t`	`float`	Parameter value in `[0, 1]`.	required
`controls_x`	`numpy.ndarray of shape (n + 1,)`	x-coordinates of the control points. Defaults to `self.controls_x`.	`None`

Returns:

Type	Description
`float`	x-coordinate of the Bézier curve at `t`.

poly_z ¶

poly_z(t: float, controls_z: ndarray = None) -> float

Evaluate the CDF of the Bézier distribution at parameter t.

Computes F(t) = B_z(t) = ∑ᵢ B_{n,i}(t) * zᵢ (Eq. 2 of the paper). Returns a value in [0, 1] when control points satisfy the CDF boundary conditions.

Parameters:

Name	Type	Description	Default
`t`	`float`	Parameter value in `[0, 1]`.	required
`controls_z`	`numpy.ndarray of shape (n + 1,)`	z-coordinates of the control points. Defaults to `self.controls_z`.	`None`

Returns:

Type	Description
`float`	CDF value at parameter `t`.

root_find ¶

root_find(x: float, method: str = 'brentq') -> float

Find the parameter t such that B_x(t) = x.

Solves poly_x(t) - x = 0 on [0, 1] using the specified root-finding method.

Parameters:

Name	Type	Description	Default
`x`	`float`	Target x-coordinate within the support `[controls_x[0], controls_x[-1]]`.	required
`method`	`(brentq, bisect)`	Root-finding algorithm. Default is `'brentq'`.	`'brentq'`

Returns:

Type	Description
`float`	Parameter `t` in `[0, 1]` satisfying `B_x(t) ≈ x`.

Raises:

Type	Description
`ValueError`	If the root-finder cannot bracket a root (`x` outside support).

eval_t ¶

eval_t(t: float) -> tuple[float, float]

Evaluate the Bezier random variable at a given parameter value t.

Parameters:

Name	Type	Description	Default
`t`	`float`	The parameter value at which to evaluate the curve (typically in [0, 1]).	required

Returns:

Type	Description
`tuple[float, float]`	A tuple containing the (x, z) coordinates of the point on the curve w.r.t. t.

eval_x ¶

eval_x(x: float) -> tuple[float, float]

Evaluate the Bezier random variable at a given x-coordinate.

Parameters:

Name	Type	Description	Default
`x`	`float`	The x-coordinate at which to evaluate the Bezier curve.	required

Returns:

Type	Description
`tuple[float, float]`	A tuple containing the (x, z) coordinates of the point on the curve w.r.t. x.

cdf_x ¶

cdf_x(x: float) -> float

Compute the cumulative distribution function (CDF) at a given x-coordinate.

Parameters:

Name	Type	Description	Default
`x`	`float`	The x-coordinate at which to evaluate the CDF.	required

Returns:

Type	Description
`float`	The CDF value at the given x-coordinate.

quantile ¶

quantile(alpha: float, method: str = 'brentq') -> float

Compute the quantile (inverse CDF) at probability level alpha.

Parameters:

Name	Type	Description	Default
`alpha`	`float`	Probability level in `[0, 1]`.	required
`method`	`(brentq, bisect)`	Root-finding algorithm used to invert the CDF. Default is `'brentq'`.	`'brentq'`

Returns:

Type	Description
`float`	`x` such that `F(x) = alpha`.

pdf_t ¶

pdf_t(t: float) -> float

Evaluate the Bézier PDF at parameter t.

Computes f(t) = B'_z(t) / B'_x(t) where the derivatives are degree-(n-1) Bézier curves in the differences Δzᵢ and Δxᵢ (Section 2 of the paper).

Parameters:

Name	Type	Description	Default
`t`	`float`	Parameter value in `[0, 1]`.	required

Returns:

Type	Description
`float`	PDF value at `t`.

pdf_x ¶

pdf_x(x: float) -> float

Compute the probability density function (PDF) of the Bezier random variable at a given x.

Parameters:

Name	Type	Description	Default
`x`	`float`	The x-coordinate at which to evaluate the PDF.	required

Returns:

Type	Description
`float`	The computed PDF value at x.

pdf_numerator_t ¶

pdf_numerator_t(t: float) -> float

Evaluate the numerator of the Bézier PDF at parameter t.

Computes B'_z(t) = ∑ᵢ B_{n-1,i}(t) * Δzᵢ, the z-derivative component of the PDF. Used internally when only the numerator is needed (e.g. for MLE gradient computations).

Parameters:

Name	Type	Description	Default
`t`	`float`	Parameter value in `[0, 1]`.	required

Returns:

Type	Description
`float`	Numerator `B'_z(t)` of the PDF at `t`.

raw_moment ¶

raw_moment(r: int) -> float

Compute and cache the raw moment E[X^r] of the distribution.

r = 1 uses the closed-form expression derived from the Bézier curve properties; r >= 2 integrates x^r * f(x) numerically over the support via :func:scipy.integrate.quad.

Parameters:

Name	Type	Description	Default
`r`	`int`	Order of the raw moment (`r >= 1`).	required

Returns:

Type	Description
`float`	`E[X^r]`.

mean ¶

mean() -> float

Mean E[X] (closed form).

variance ¶

variance() -> float

Variance E[X^2] - E[X]^2.

skewness ¶

skewness() -> float

Skewness E[(X - μ)^3] / σ^3.

kurtosis ¶

kurtosis() -> float

Kurtosis E[(X - μ)^4] / σ^4 (Pearson, not excess).

random ¶

random(n_sims: int, *, rng: Generator | int | None = None)

Generate random samples from the Bezier random variable.

This method generates n_sims random samples from the Bezier random variable by evaluating the inverse CDF at uniformly distributed random values in the interval [0, 1].

Parameters:

Name	Type	Description	Default
`n_sims`	`int`	The number of random samples to generate.	required
`rng`	`Generator \| int \| None`	Pseudorandom-number generator state. If None (default), a new `numpy.random.Generator` is created with fresh OS entropy. Any value accepted by :func:`numpy.random.default_rng` (e.g. a seed integer or :class:`~numpy.random.SeedSequence`) is also allowed.	`None`

Returns:

Type	Description
`(ndarray, shape(n_sims))`	Random samples drawn from the Bézier distribution.

plot_cdf ¶

plot_cdf(data: ndarray = None, num_points: int = 100, ax: Axes = None, show: bool = True)

Plot the Bézier CDF, optionally alongside the empirical CDF.

Parameters:

Name	Type	Description	Default
`data`	`array - like`	Observed sample. If provided, the empirical CDF is overlaid. If `None`, a linspace over the support is used.	`None`
`num_points`	`int`	Number of evaluation points when `data` is `None`. Default is `100`.	`100`
`ax`	`Axes`	Axes to plot on. Defaults to `plt.gca()`.	`None`
`show`	`bool`	If `True` (default), call `plt.show()` after plotting.	`True`

plot_pdf ¶

plot_pdf(data: ndarray = None, num_points: int = 100, ax: Axes = None, show: bool = True)

Plot the Bézier PDF.

Parameters:

Name	Type	Description	Default
`data`	`array - like`	Points at which to evaluate the PDF. If `None`, a linspace over the support is used.	`None`
`num_points`	`int`	Number of evaluation points when `data` is `None`. Default is `100`.	`100`
`ax`	`Axes`	Axes to plot on. Defaults to `plt.gca()`.	`None`
`show`	`bool`	If `True` (default), call `plt.show()` after plotting.	`True`

Quick Example

from bezierv.classes.bezierv import Bezierv
import numpy as np

# Create a simple Bézier distribution
controls_x = [0, 1, 2, 3]
controls_z = [0, 0.3, 0.8, 1]

bezier_rv = Bezierv(3, controls_x, controls_z)

# Evaluate PDF and CDF
pdf_value = bezier_rv.pdf_x(1.5)
cdf_value = bezier_rv.cdf_x(1.5)

# Generate random samples
samples = bezier_rv.random(100, rng=42)

InteractiveBezierv - Interactive Editor¶

Interactive Bokeh-based editor for hands-on exploration of Bézier distributions.

bezierv.classes.bezierv.InteractiveBezierv ¶

InteractiveBezierv(controls_x, controls_z)

Manages an interactive Bezier distribution in a Bokeh plot.

Interactive Session

from bezierv.classes.bezierv import InteractiveBezierv
from bokeh.plotting import curdoc

# Start with uniform distribution
controls_x = [0, 1, 2, 3]
controls_z = [0, 0.33, 0.67, 1]

editor = InteractiveBezierv(controls_x, controls_z)
curdoc().add_root(editor.layout)

# Run with: bokeh serve --show script.py

Distribution Fitting¶

DistFit - Fit Distributions to Data¶

Fits Bézier distributions to empirical data using various optimization algorithms.

bezierv.classes.distfit.DistFit ¶

DistFit(data: List, n: int = 10, init_x: array = None, init_z: array = None, init_t: array = None, init_w: array = None, emp_cdf_data: array = None, method_init_x: str = 'quantile')

A class to fit a Bezier random variable to empirical data.

Attributes:

Name	Type	Description
`data`	`List`	The empirical data to fit the Bezier random variable to.
`n`	`int`	The number of control points minus one for the Bezier curve.
`init_x`	`array`	Initial control points for the x-coordinates of the Bezier curve.
`init_z`	`array`	Initial control points for the z-coordinates of the Bezier curve.
`init_t`	`array`	Initial parameter values.
`init_w`	`array`	Initial weight vector on the probability simplex.
`emp_cdf_data`	`array`	The empirical cumulative distribution function (CDF) data derived from the empirical data.
`bezierv`	`Bezierv`	An instance of the Bezierv class representing the Bezier random variable.
`m`	`int`	The number of empirical data points.
`mse`	`float`	The mean squared error of the fit, initialized to infinity.
`nll`	`float`	The negative log-likelihood of the fit, initialized to infinity.

Initialize the DistFit class with empirical data and parameters for fitting a Bezier random variable.

Parameters:

Name	Type	Description	Default
`data`	`List`	The empirical data to fit the Bezier random variable to.	required
`n`	`int`	The number of control points minus one for the Bezier curve (default is 10).	`10`
`init_x`	`array`	Initial control points for the x-coordinates of the Bezier curve (default is None).	`None`
`init_z`	`array`	Initial control points for the z-coordinates of the Bezier curve (default is None).	`None`
`init_t`	`array`	Initial parameter values (default is None).	`None`
`init_w`	`array`	Initial weight vector on the probability simplex (default is None).	`None`
`emp_cdf_data`	`array`	The empirical cumulative distribution function (CDF) data derived from the empirical data (default is None).	`None`
`method_init_x`	`str`	Method to initialize the x-coordinates of the control points (default is 'quantile').	`'quantile'`

Methods:¶

fit ¶

fit(method: str = 'mse', algorithm: str = 'projected_gradient', options: FitOptions = None) -> tuple[Bezierv, float]

Fit the bezierv distribution to the data.

Parameters:

Name	Type	Description	Default
`method`	`str`	The fitting method to use. Options are 'mse' for mean squared error or 'mle' for maximum likelihood estimation (default is 'mse').	`'mse'`
`algorithm`	`str`	The fitting algorithm to use when `method='mse'`. Options are 'projected_gradient', 'solver', or 'nelder_mead'. Ignored when `method='mle'`. Default is 'projected_gradient'. The legacy names 'projgrad', 'nonlinear', and 'neldermead' are still accepted as deprecated aliases.	`'projected_gradient'`
`options`	`FitOptions`	Algorithm-specific options. Must match the chosen method/algorithm: `ProjGradOptions` for ('mse', 'projected_gradient'), `NonLinearOptions` for ('mse', 'solver'), `NelderMeadOptions` for ('mse', 'nelder_mead'), `MLEOptions` for ('mle'). Defaults to the matching dataclass with its default field values.	`None`

Returns:

Name	Type	Description
`bezierv`	`Bezierv`	A copy of the fitted Bezierv instance with updated control points.
`metric`	`float`	Final value of the objective: mean squared error when `method='mse'`, or negative log-likelihood when `method='mle'`.

get_controls_z ¶

get_controls_z() -> array

Compute the control points for the z-coordinates of the Bezier curve.

Returns:

Type	Description
`array`	The control points for the z-coordinates of the Bezier curve, evenly spaced between 0 and 1.

get_controls_x ¶

get_controls_x(method: str) -> array

Compute the control points for the x-coordinates of the Bezier curve.

'quantile' method is used to determine the control points based on the data quantiles.

Parameters:

Name	Type	Description	Default
`method`	`str`	The method to use for initializing the x-coordinates of the control points.	required

Returns:

Type	Description
`array`	The control points for the x-coordinates of the Bezier curve.

Fitting Example

from bezierv.classes.distfit import DistFit, ProjGradOptions
import numpy as np

# Generate sample data
data = np.random.gamma(2, 3, 500)

# Fit Bézier distribution (defaults to MSE + projected gradient)
# The step size defaults to η = 1/β with β = (2/m) λ_max(BᵀB), the
# Lipschitz constant of the MSE gradient — no manual tuning needed.
fitter = DistFit(data, n=5)
bezier_rv, mse = fitter.fit(algorithm="projected_gradient")

# Or override the automatic step size and other tunables
bezier_rv, mse = fitter.fit(
    algorithm="projected_gradient",
    options=ProjGradOptions(step_size=5e-4, max_iter=2000, threshold=1e-4),
)

print(f"Fit quality (MSE): {mse:.6f}")

# Visualize results
bezier_rv.plot_cdf(data)
bezier_rv.plot_pdf()

Fit Options¶

Algorithm-specific tunables for DistFit.fit are grouped into dataclasses. Pass an instance of the matching dataclass as the options argument; omit it to use the defaults.

Method + Algorithm	Options class	Fields (defaults)
`mse` + `projected_gradient`	`ProjGradOptions`	`step_size=None` (auto: `η = 1/β`, `β = (2/m) λ_max(BᵀB)`), `max_iter=1000`, `threshold=1e-3`
`mse` + `solver`	`NonLinearOptions`	`solver='ipopt'`, `solver_options={'timelimit': 60, 'tee': False}`
`mse` + `nelder_mead`	`NelderMeadOptions`	`max_iter=1000`
`mle`	`MLEOptions`	`max_iter=1000`, `tol=1e-3`, `tol_res_root=1e-5`, `tol_lambda_root=1e-5`, `max_iters_root=100`

from bezierv.classes.distfit import (
    DistFit,
    ProjGradOptions, NonLinearOptions, NelderMeadOptions, MLEOptions,
)

fitter = DistFit(data, n=5)

# Defaults
fitter.fit()

# Nonlinear with IPOPT options; solver_options is forwarded as
# ``pyo_solver.solve(model, **solver_options)`` — Pyomo ``solve`` kwargs
# (``timelimit``, ``tee``) sit at the top level, solver-native options
# under ``options``.
fitter.fit(
    algorithm="solver",
    options=NonLinearOptions(
        solver="ipopt",
        solver_options={"timelimit": 120, "tee": False, "options": {"max_iter": 5000, "tol": 1e-8}},
    ),
)

# MLE with tighter tolerances
fitter.fit(method="mle", options=MLEOptions(tol=1e-6, max_iter=5000))

Options type must match the algorithm

Passing ProjGradOptions while selecting algorithm="solver" raises TypeError. The options class encodes the algorithm choice.

Convolution Operations¶

Convolver - Sum of Random Variables¶

Computes convolutions (sums) of independent Bézier random variables using Monte Carlo simulation or exact numerical integration.

bezierv.classes.convolver.Convolver ¶

Convolver(list_bezierv: List[Bezierv])

Convolution of independent Bézier random variables.

Supports Monte Carlo convolution for any number of variables and exact convolution via numerical integration for the two-variable case.

Attributes:

Name	Type	Description
`list_bezierv`	`List[Bezierv]`	Bézier random variables to be convolved. Each element is validated (lengths, ordering, initialization) at construction time.

Initialize a Convolver from a list of Bézier random variables.

Parameters:

Name	Type	Description	Default
`list_bezierv`	`List[Bezierv]`	Bézier random variables to be convolved.	required

Raises:

Type	Description
`ValueError`	If any element fails the length, ordering, or initialization checks performed by :class:`~bezierv.classes.bezierv.Bezierv`.

Methods:¶

convolve ¶

convolve(n_sims: int = 1000, *, rng: Generator | int | None = None, **kwargs: Any) -> tuple[Bezierv, float]

Convolve the Bezier RVs via Monte Carlo and fit a Bezierv to the sum.

Draws n_sims samples from each Bézier RV in :attr:list_bezierv using a shared PRNG stream, sums them, and fits a new Bézier random
variable to the resulting sample via :class:DistFit.

Parameters:

Name	Type	Description	Default
`n_sims`	`int`	Number of Monte Carlo samples. Default is `1000`.	`1000`
`rng`	`Generator \| int \| None`	Shared PRNG stream for all sampling.	`None`
`**kwargs`	`Any`	Init options for DistFit(...): n, init_x, init_z, init_t, emp_cdf_data, method_init_x Fit options for DistFit.fit(...): method, algorithm, options	`{}`

Returns:

Name	Type	Description
`bezierv`	`Bezierv`	Bézier random variable fitted to the convolved sample.
`metric`	`float`	Final objective value reported by :meth:`DistFit.fit` (MSE for `method='mse'`, NLL for `method='mle'`).

Raises:

Type	Description
`TypeError`	If `kwargs` contains a key that is not recognized as either a `DistFit` init option or a `DistFit.fit` option.

convolve_exact ¶

convolve_exact(n_points: int = 1000, **kwargs) -> tuple[Bezierv, float]

Perform convolution of two Bezier random variables using numerical integration.

This method generates a grid of points and computes the exact CDF at each point using the exact_cdf_two_bezierv method, then fits a Bezierv to the resulting CDF values.

Parameters:

Name	Type	Description	Default
`n_points`	`int`	Number of points to evaluate the exact CDF (default is 1000).	`1000`
`**kwargs`		Init options for DistFit(...): n, init_x, init_z, init_t, emp_cdf_data, method_init_x Fit options for DistFit.fit(...): method, algorithm, options	`{}`

Returns:

Name	Type	Description
`bezierv`	`Bezierv`	The fitted Bezierv representing the convolution.
`metric`	`float`	Final objective value reported by :meth:`DistFit.fit` (MSE for `method='mse'`, NLL for `method='mle'`).

Raises:

Type	Description
`ValueError`	If the number of Bezier RVs is not exactly 2.

exact_cdf_two_bezierv ¶

exact_cdf_two_bezierv(z: float) -> float

Compute the exact CDF of the convolution Z = X + Y at point z using numerical integration.

Parameters:

Name	Type	Description	Default
`z`	`float`	The point at which to evaluate the CDF.	required

Returns:

Type	Description
`float`	The CDF value F_Z(z) at point z.

Method Comparison

We recommend using monte carlo convolution.

Convolution Example

from bezierv.classes.convolver import Convolver
from bezierv.classes.distfit import DistFit
import numpy as np

# Create two distributions
data1 = np.random.gamma(2, 2, 1000)  # Task 1 duration
data2 = np.random.lognormal(1, 0.5, 1000)  # Task 2 duration

# Fit distributions
rv1, _ = DistFit(data1, n=5).fit()
rv2, _ = DistFit(data2, n=5).fit()

# Compute convolution
convolver = Convolver([rv1, rv2])

total_mc, _ = convolver.convolve(n_sims=1000, rng=42)

print(f"Expected total time (MC): {total_mc.mean():.2f}")

Optimization Algorithms¶

These algorithms are used internally by DistFit to optimize the fitting process. You typically don't need to call these directly.

Projected Gradient Method¶

Fast and stable gradient-based optimization with projection onto the feasible region.

bezierv.algorithms.proj_grad ¶

Classes¶

Functions:¶

compute_bernstein_basis ¶

compute_bernstein_basis(n: int, t: ndarray, bezierv: Bezierv) -> ndarray

Precompute the degree-n Bernstein basis matrix for the Bézier CDF.

Evaluates all n + 1 degree-n Bernstein basis polynomials at each parameter value in t. The resulting matrix B satisfies B @ z = F̂ in the minimum-error objective (Section 3 of the paper). Unlike the basis in primal_grad, this uses degree n (for the CDF) rather than degree n - 1 (for the PDF).

Parameters:

Name	Type	Description	Default
`n`	`int`	Degree of the Bézier curve (`n + 1` control points).	required
`t`	`(ndarray, shape(m))`	Parameter values in `[0, 1]` at which to evaluate the basis.	required
`bezierv`	`Bezierv`	Bézier random variable providing the precomputed binomial coefficients `comb`.	required

Returns:

Type	Description
`(ndarray, shape(m, n + 1))`	Basis matrix where `B[j, i]` is the `i`-th degree-`n` Bernstein basis polynomial evaluated at `t[j]`.

lipschitz_mse ¶

lipschitz_mse(B: ndarray) -> float

Compute the Lipschitz constant of the MSE gradient for projected gradient descent.

Returns β = (2/m) * λ_max(BᵀB), the smallest Lipschitz constant of ∇f for the normalized MSE objective f(z) = (1/m) ‖Bz - F̂‖² (Section 3 of the paper), where B is the Bernstein basis matrix from :func:compute_bernstein_basis. The projected gradient step η = 1/β is the largest constant step that still guarantees monotone descent on this smooth convex objective.

Two branches are selected by the number of control points n + 1 = B.shape[1]: for n + 1 <= 256 the exact value is computed from σ_max(B) = np.linalg.norm(B, 2) so that BᵀB is never formed (avoiding the condition-number squaring of an eigvalsh call); for n + 1 > 256 the dominant eigenvalue is estimated by matrix-free Lanczos (:func:scipy.sparse.linalg.eigsh on a :class:~scipy.sparse.linalg.LinearOperator implementing the matvec v -> Bᵀ(Bv)), then inflated by a factor of 1.01 because Lanczos Ritz values lower-bound λ_max, so the inflation ensures β over-bounds the true Lipschitz constant and η = 1/β stays safe.

Parameters:

Name	Type	Description	Default
`B`	`(ndarray, shape(m, n + 1))`	Bernstein basis matrix. Each row holds the `n + 1` degree-`n` Bernstein basis polynomials evaluated at one parameter `t[j]`.	required

Returns:

Type	Description
`float`	`β = (2/m) * λ_max(BᵀB)`, safety-inflated on the Lanczos branch.

Raises:

Type	Description
`ValueError`	If `B` has zero rows (`m = 0`) or the largest singular value of `B` is zero (degenerate basis).

grad ¶

grad(n: int, t: ndarray, bezierv: Bezierv, controls_z: ndarray, emp_cdf_data: ndarray, basis_matrix: ndarray = None) -> ndarray

Compute the gradient of the MSE objective w.r.t. the z control points.

Evaluates ∇_z (1/m) ‖Bz - F̂‖² = (2/m) Bᵀ(Bz - F̂) where B is the Bernstein basis matrix and F̂ are the empirical CDF values (Section 3 of the paper). The 1/m normalization matches the MSE objective and is the gradient whose Lipschitz constant is returned by :func:lipschitz_mse.

Parameters:

Name	Type	Description	Default
`n`	`int`	Degree of the Bézier curve.	required
`t`	`(ndarray, shape(m))`	Parameter values in `[0, 1]` for each observation.	required
`bezierv`	`Bezierv`	Bézier random variable providing the Bernstein basis evaluator.	required
`controls_z`	`(ndarray, shape(n + 1))`	Current z-coordinates of the control points.	required
`emp_cdf_data`	`(ndarray, shape(m))`	Empirical CDF values at each observation.	required
`basis_matrix`	`(ndarray, shape(m, n + 1))`	Precomputed basis matrix from :func:`compute_bernstein_basis`. Computed on the fly if not provided.	`None`

Returns:

Type	Description
`(ndarray, shape(n + 1))`	Gradient `(2/m) Bᵀ(Bz - F̂)`.

project_z ¶

project_z(controls_z: ndarray) -> ndarray

Project the z control points onto the bounded monotone feasible set.

Delegates to :func:~bezierv.algorithms.isotonic_reg.project with boundary conditions z₀ = 0 and zₙ = 1 (CDF constraints from Section 3 of the paper).

Parameters:

Name	Type	Description	Default
`controls_z`	`(ndarray, shape(n + 1))`	Current z-coordinates of the control points.	required

Returns:

Type	Description
`(ndarray, shape(n + 1))`	Projected control points satisfying `0 = z[0] <= ... <= z[n] = 1`.

fit ¶

fit(n: int, bezierv: Bezierv, init_x: ndarray, init_z: ndarray, t: ndarray, emp_cdf_data: ndarray, step_size: float = None, max_iter: int = 1000, threshold: float = 0.001) -> tuple[Bezierv, float]

Fit a Bézier distribution via projected gradient descent on the MSE objective.

Minimizes (1/m) ‖Bz - F̂‖² over the bounded monotone feasible set by iterating a gradient step followed by projection via :func:project_z. Implements Algorithm 1 from Section 3 of the paper.

Parameters:

Name	Type	Description	Default
`n`	`int`	Degree of the Bézier curve (`n + 1` control points).	required
`bezierv`	`Bezierv`	Bézier random variable to be updated with the fitted parameters.	required
`init_x`	`(ndarray, shape(n + 1))`	Fixed x-coordinates of the control points (support knots).	required
`init_z`	`(ndarray, shape(n + 1))`	Initial z-coordinates of the control points.	required
`t`	`(ndarray, shape(m))`	Pre-computed parameter values for each observation, from :func:`~bezierv.algorithms.utils.get_t`.	required
`emp_cdf_data`	`(ndarray, shape(m))`	Empirical CDF values at each observation.	required
`step_size`	`float`	Gradient descent step size `η`. If `None` (default), an automatic step `η = 1/β` is used, where `β = (2/m) λ_max(BᵀB)` is computed by :func:`lipschitz_mse`. Pass a positive float to override the automatic choice.	`None`
`max_iter`	`int`	Maximum number of projected gradient iterations. Default `1000`.	`1000`
`threshold`	`float`	Convergence tolerance on the relative change `‖z_{k+1} - z_k‖ / max(1, ‖z_k‖)`. Default `1e-3`.	`0.001`

Returns:

Name	Type	Description
`bezierv`	`Bezierv`	Updated Bézier random variable with fitted control points.
`mse`	`float`	Final mean squared error between the fitted CDF and the empirical CDF.

When to Use

Default choice for most fitting problems. Provides good balance of speed and accuracy.

Nonlinear Optimization¶

Pyomo-based nonlinear program solved with an external NLP solver (IPOPT by default). Recommended when fit accuracy matters more than speed and a suitable solver is available on the system.

bezierv.algorithms.non_linear ¶

Classes¶

Functions:¶

fit ¶

fit(n: int, m: int, data: ndarray, bezierv: Bezierv, init_x: ndarray, init_z: ndarray, init_t: ndarray, emp_cdf_data: ndarray, solver: str, solver_options: dict = None, feas_tol: float = 0.0001) -> tuple[Bezierv, float]

Fit a Bézier distribution via an external solver (e.g. IPOPT).

Formulates the minimum-error estimation problem as a Pyomo ConcreteModel and solves it with the specified NLP solver.

Parameters:

Name	Type	Description	Default
`n`	`int`	Degree of the Bézier curve (`n + 1` control points).	required
`m`	`int`	Sample size.	required
`data`	`(ndarray, shape(m))`	Sorted observed sample values.	required
`bezierv`	`Bezierv`	Bézier random variable to be updated with the fitted parameters.	required
`init_x`	`(ndarray, shape(n + 1))`	Initial x-coordinates of the control points.	required
`init_z`	`(ndarray, shape(n + 1))`	Initial z-coordinates (CDF values) of the control points.	required
`init_t`	`(ndarray, shape(m))`	Initial parameter values `t ∈ [0, 1]` for each observation.	required
`emp_cdf_data`	`(ndarray, shape(m))`	Empirical CDF values at each observation.	required
`solver`	`str`	Name of the NLP solver recognized by Pyomo (e.g. `'ipopt'`).	required
`solver_options`	`dict`	Keyword arguments forwarded to `pyo_solver.solve(model, **solver_options)`. Supports Pyomo's `solve` kwargs (e.g. `tee`, `timelimit`) and solver-native options nested under the `options` key (e.g. `{'timelimit': 60, 'tee': False, 'options': {'max_iter': 5000, 'tol': 1e-8}}` for IPOPT). Defaults to `None`.	`None`
`feas_tol`	`float`	Maximum allowed constraint or bound violation when accepting an incumbent solution from a non-optimal termination (e.g. `maxIterations`, `maxTimeLimit`). If the worst violation exceeds this tolerance, a :class:`RuntimeError` is raised. Defaults to `1e-4`, matching IPOPT's `constr_viol_tol` default.	`0.0001`

Returns:

Name	Type	Description
`bezierv`	`Bezierv`	Updated Bézier random variable. Unchanged if the solver does not reach an optimal solution.
`mse`	`float`	Final MSE at the optimal solution, or `nan` if the solver fails.

Raises:

Type	Description
`RuntimeError`	If the solver returns a non-optimal acceptable termination (`maxIterations`, `maxTimeLimit`, `locallyOptimal`, `feasible`) but the incumbent violates the constraints or variable bounds by more than `feas_tol`, or if the solver terminates with any other (non-acceptable) status.
`Exception`	Any exception raised by the Pyomo solver call is propagated to the caller unchanged.

Notes

Constraints enforce: boundary conditions z₀ = 0, zₙ = 1; monotonicity of x and z; support bounds x₀ ≤ data[0], xₙ ≥ data[-1]; and parameter boundary t₁ = 0, tₘ = 1.

When to Use

Use when accuracy is more important than runtime, you can afford failures, and a Pyomo-compatible NLP solver such as IPOPT is installed.

Requires an external solver

The default solver is ipopt. You must have it installed and available on your PATH (or pass NonLinearOptions(solver=...) to use a different Pyomo-compatible NLP solver).

Nelder-Mead¶

Derivative-free optimization method that doesn't require gradients.

bezierv.algorithms.nelder_mead ¶

Classes¶

Functions:¶

objective_function ¶

objective_function(concatenated: ndarray, n: int, m: int, data: ndarray, bezierv: Bezierv, emp_cdf_data: ndarray) -> float

Compute the MSE between the fitted Bézier CDF and the empirical CDF.

Evaluates the minimum-error objective (1/m) ∑ⱼ (F(xⱼ) - F̂ⱼ)² where F is the Bézier CDF parameterized by the given control points and F̂ⱼ are the empirical CDF values (Section 3 of the paper).

Parameters:

Name	Type	Description	Default
`concatenated`	`(ndarray, shape(2 * (n + 1)))`	Concatenated control point vector `[x₀, …, xₙ, z₀, …, zₙ]`.	required
`n`	`int`	Degree of the Bézier curve (`n + 1` control points).	required
`m`	`int`	Sample size.	required
`data`	`(ndarray, shape(m))`	Sorted observed sample values.	required
`bezierv`	`Bezierv`	Bézier random variable providing the CDF evaluator `poly_z`.	required
`emp_cdf_data`	`(ndarray, shape(m))`	Empirical CDF values at each observation.	required

Returns:

Type	Description
`float`	Mean squared error between the fitted CDF and the empirical CDF.

objective_function_lagrangian ¶

objective_function_lagrangian(concatenated: ndarray, n: int, m: int, data: ndarray, bezierv: Bezierv, emp_cdf_data: ndarray, penalty_weight: float = 1000.0) -> float

Compute the penalized MSE objective for unconstrained Nelder-Mead optimization.

Augments :func:objective_function with a Lagrangian penalty that enforces the Bézier feasibility constraints — boundary conditions z₀ = 0, zₙ = 1, monotonicity of z and x, and support matching x₀ = data[0], xₙ = data[-1] — since Nelder-Mead does not support constraints natively (Section 3 of the paper, Wagner baseline method).

Parameters:

Name	Type	Description	Default
`concatenated`	`(ndarray, shape(2 * (n + 1)))`	Concatenated control point vector `[x₀, …, xₙ, z₀, …, zₙ]`.	required
`n`	`int`	Degree of the Bézier curve (`n + 1` control points).	required
`m`	`int`	Sample size.	required
`data`	`(ndarray, shape(m))`	Sorted observed sample values.	required
`bezierv`	`Bezierv`	Bézier random variable providing the CDF evaluator `poly_z`.	required
`emp_cdf_data`	`(ndarray, shape(m))`	Empirical CDF values at each observation.	required
`penalty_weight`	`float`	Scaling factor for the constraint violation penalty. Default is `1e3`.	`1000.0`

Returns:

Type	Description
`float`	Penalized objective `MSE + penalty_weight * constraint_violation`. Returns `inf` if root-finding fails for the given `x` knots.

fit ¶

fit(n: int, m: int, data: ndarray, bezierv: Bezierv, init_x: ndarray, init_z: ndarray, emp_cdf_data: ndarray, max_iter: int) -> tuple[Bezierv, float]

Fit a Bézier distribution using the Nelder-Mead derivative-free method.

Minimizes :func:objective_function_lagrangian via scipy.optimize.minimize with method='Nelder-Mead'. This is the derivative-free baseline (Wagner's PRIME approach) benchmarked against the first-order methods in Section 7 of the paper.

Parameters:

Name	Type	Description	Default
`n`	`int`	Degree of the Bézier curve (`n + 1` control points).	required
`m`	`int`	Sample size.	required
`data`	`(ndarray, shape(m))`	Sorted observed sample values.	required
`bezierv`	`Bezierv`	Bézier random variable to be updated with the fitted parameters.	required
`init_x`	`(ndarray, shape(n + 1))`	Initial x-coordinates of the control points.	required
`init_z`	`(ndarray, shape(n + 1))`	Initial z-coordinates (CDF values) of the control points.	required
`emp_cdf_data`	`(ndarray, shape(m))`	Empirical CDF values at each observation.	required
`max_iter`	`int`	Maximum number of Nelder-Mead iterations.	required

Returns:

Name	Type	Description
`bezierv`	`Bezierv`	Updated Bézier random variable with fitted control points.
`mse`	`float`	Final mean squared error of the fit, evaluated via :func:`objective_function`.

Algorithm Utilities¶

Helper functions and utilities used by the optimization algorithms.

bezierv.algorithms.utils ¶

Classes¶

Functions:¶

root_find ¶

root_find(n: int, bezierv: Bezierv, controls_x: ndarray, data_point: float) -> float

Find the parameter t such that the Bézier x-coordinate equals a given data point.

Uses Brent's method to solve B_x(t) - data_point = 0 on [0, 1], where B_x(t) is the x-component of the degree-n Bézier curve (Eq. 2 of the paper).

Parameters:

Name	Type	Description	Default
`n`	`int`	Degree of the Bézier curve. A degree-`n` curve has `n + 1` control points.	required
`bezierv`	`Bezierv`	Fitted Bézier random variable providing the Bernstein basis evaluator.	required
`controls_x`	`(ndarray, shape(n + 1))`	x-coordinates of the `n + 1` control points.	required
`data_point`	`float`	Observed value for which the inverse parameter `t` is sought.	required

Returns:

Type	Description
`float`	Parameter value `t` in `[0, 1]` satisfying `B_x(t) = data_point`.

Raises:

Type	Description
`ValueError`	If `brentq` cannot bracket a root, i.e. `data_point` lies outside the range of `B_x` on `[0, 1]`.

get_t ¶

get_t(n: int, m: int, data: ndarray, bezierv: Bezierv, controls_x: ndarray) -> ndarray

Compute the Bézier parameter t for each data point via root-finding.

For each observation in data, solves B_x(t) = data[i] using :func:root_find to obtain the corresponding parameter value on [0, 1].

Parameters:

Name	Type	Description	Default
`n`	`int`	Degree of the Bézier curve (`n + 1` control points).	required
`m`	`int`	Number of observations (`len(data)`).	required
`data`	`(ndarray, shape(m))`	Observed sample values in `[a, b]`.	required
`bezierv`	`Bezierv`	Fitted Bézier random variable providing the Bernstein basis evaluator.	required
`controls_x`	`(ndarray, shape(n + 1))`	x-coordinates of the `n + 1` control points.	required

Returns:

Type	Description
`(ndarray, shape(m))`	Parameter values `t[i]` in `[0, 1]` such that `B_x(t[i]) = data[i]` for each `i`.

API Reference¶

Core Classes¶

Bezierv - Main Distribution Class¶

bezierv.classes.bezierv.Bezierv ¶

Methods:¶

update_bezierv ¶

combinations ¶

deltas ¶

bernstein ¶

poly_x ¶

poly_z ¶

root_find ¶

eval_t ¶

eval_x ¶

cdf_x ¶

quantile ¶

pdf_t ¶

pdf_x ¶

pdf_numerator_t ¶

raw_moment ¶

mean ¶

variance ¶

skewness ¶

kurtosis ¶

random ¶

plot_cdf ¶

plot_pdf ¶

InteractiveBezierv - Interactive Editor¶

bezierv.classes.bezierv.InteractiveBezierv ¶

Distribution Fitting¶

DistFit - Fit Distributions to Data¶

bezierv.classes.distfit.DistFit ¶

Methods:¶

fit ¶

get_controls_z ¶

get_controls_x ¶

Fit Options¶

Convolution Operations¶

Convolver - Sum of Random Variables¶

bezierv.classes.convolver.Convolver ¶

Methods:¶

convolve ¶

convolve_exact ¶

exact_cdf_two_bezierv ¶

Optimization Algorithms¶

Projected Gradient Method¶

bezierv.algorithms.proj_grad ¶

Classes¶

Functions:¶

compute_bernstein_basis ¶

lipschitz_mse ¶

grad ¶

project_z ¶

fit ¶

Nonlinear Optimization¶

bezierv.algorithms.non_linear ¶

Classes¶

Functions:¶

fit ¶

Nelder-Mead¶

bezierv.algorithms.nelder_mead ¶

Classes¶

Functions:¶

objective_function ¶

objective_function_lagrangian ¶

fit ¶

Algorithm Utilities¶

bezierv.algorithms.utils ¶

Classes¶

Functions:¶

root_find ¶

get_t ¶

See Also¶