freealg.AlgebraicForm.fit#

AlgebraicForm.fit(deg_m, deg_z, reg=0.0, r_min=1.8, r_max=2.2, n_r=5, y_scale=1.0, gamma=1.0, n_samples=4096, cut_eps=0.01, triangular=None, mu='auto', mu_reg=None, normalize=True, verbose=False, plot=False)#

Fit an algebraic structure to the input data.

Parameters:

deg_mint

\(\mathrm{deg}_m(P)\): degree polynomial \(P(z, m)\) in \(m\).

deg_zint

\(\mathrm{deg}_z(P)\): degree polynomial \(P(z, m)\) in \(z\).

regfloat, default=0.0

Tikhonov regularization parameter for fitting.

r_minarray_like, default=1.8

Minimum radius in Joukowski plane where contours around support intervals are generated.

r_maxarray_like, default=1.8

Maximum radius in Joukowski plane where contours around support intervals are generated.

n_rarray_like, default=1.8

Number of radii in Joukowski plane where contours around support intervals are generated.

y_scalefloat, default=1.0

The scale of y for sampling points in z plane.

gammafloat, default=1.0

The exponent in which y components of sampling points increase ( only in log-scale).

n_samplesint, default=4.96

Number of angular samples on each Bernstein ellipse from \(\theta = 0\) to \(\theta = 2 \pi\).

cut_epsfloat, default=0.01

Fraction of cut width used as the distance threshold.

log=False: threshold is cut_eps * (b - a), where (b - a) is the linear width of the cut interval.
log=True: threshold is cut_eps * log(b / a), where log(b / a) is the logarithmic width of the cut interval.

triangular{None, tuple}, default=None

Structure of the coefficient index set.

None: full rectangular index set.
(a, b): banded support defined by \(a \le j - i \le b\), where each of a and b may be an integer or None.
- a=None means there is no lower bound on \(j-i\).
- b=None means there is no upper bound on \(j-i\).

Examples:

(0, None) keeps terms with \(j \ge i\).
(-E, None) reproduces the old upper-Hessenberg case \(i \le j + E\).
(a, b) with both finite keeps only the two-sided band \(a \le j-i \le b\).

muarray_like, default= 'auto'

Constraint to fit polynomial coefficients based on moments:

If an array \([\mu_0, \mu_1\), dots, mu_r]` is given, it enforces the first \(r+1\) moments. Note that \(\mu_0 = 1\) to ensure unit mass.
If instead this option is set to 'auto', and the input A is a matrix, it automatically uses the first two moments of the eigenvalues of the input matrix as moment constraints.
If None, no constraint is used.

See also

support
branch_points
atoms
stieltjes

Notes

The Stieltjes transform \(m = m(z)\) is the root of the implicit relation \(P(z, m) = 0\), where \(P\) is given by

\[P(z, m) = \sum_{i=1}^{d_z} \sum_{j=1}^{d_m} c_{ij} z^i m^j,\]

and \(d_z\), \(d_m\) are the degrees of \(P\) in \(z\) and \(m\), respectively.

This function returns the \((d_z \times d_m)\) matrix \(\mathbf{C} = [c_{ij}]\)

\[\begin{split}\mathbf{C} = \begin{bmatrix} c_{11} & \dots & c_{1 d_m} \\ \vdots & & \vdots \\ c_{d_z 1} & \dots & c_{d_z d_m} \end{bmatrix}.\end{split}\]

These coefficients are real, however, they are returned as complex numbers with their imaginary part all being zero.

This matrix is empirically fitted based on the samples of the Stieltjes transform. To get this matrix, the fit() function should have been called first.

Examples

>>> # Create a distribution
>>> from freealg.distributions import FreeLevy

>>> # Create an object of the class
>>> fl = FreeLevy(t=[2.0, 5.5], w=[0.75, 1-0.75], lam=0.1, a=0,
...               sigma=0.9)

>>> # Create an algebraic form object from the distribution
>>> from freealg import AlgebraicForm
>>> af = AlgebraicForm(fl)

>>> # Fit polynomial
>>> af.fit(deg_m=4, deg_z=1)

>>> # Get the empirically fitted polynomial coefficients
>>> print(af.coeffs.real)
[[ 1.      7.2125 12.15   16.875  24.75  ]
 [ 0.      1.      7.5    11.      0.    ]]

>>> # Compare with the exact polynomial from the distribution law
>>> print(fl.poly().real)
[[ 1.      7.2125 12.15   16.875  24.75  ]
 [ 0.      1.      7.5    11.     -0.    ]]