plementation does nothing to prevent this.

For either method,
the returned answer is not guaranteed to be globally optimal; it
may only be locally optimal, or the optimization may fail altogether.
If the data contain any of ``np.nan``, ``np.inf``, or ``-np.inf``,
the `fit` method will raise a ``RuntimeError``.

When passing a ``CensoredData`` instance to ``data``, the log-likelihood
function is defined as:

.. math::

    l(\pmb{\theta}; k) & = \sum
                            \log(f(k_u; \pmb{\theta}))
                        + \sum
                            \log(F(k_l; \pmb{\theta})) \\
                        & + \sum
                            \log(1 - F(k_r; \pmb{\theta})) \\
                        & + \sum
                            \log(F(k_{\text{high}, i}; \pmb{\theta})
                            - F(k_{\text{low}, i}; \pmb{\theta}))

where :math:`f` and :math:`F` are the pdf and cdf, respectively, of the
function being fitted, :math:`\pmb{\theta}` is the parameter vector,
:math:`u` are the indices of uncensored observations,
:math:`l` are the indices of left-censored observations,
:math:`r` are the indices of right-censored observations,
subscripts "low"/"high" denote endpoints of interval-censored observations, and
:math:`i` are the indices of interval-censored observations.

Examples
--------

Generate some data to fit: draw random variates from the `beta`
distribution

>>> import numpy as np
>>> from scipy.stats import beta
>>> a, b = 1., 2.
>>> rng = np.random.default_rng(172786373191770012695001057628748821561)
>>> x = beta.rvs(a, b, size=1000, random_state=rng)

Now we can fit all four parameters (``a``, ``b``, ``loc`` and
``scale``):

>>> a1, b1, loc1, scale1 = beta.fit(x)
>>> a1, b1, loc1, scale1
(1.0198945204435628, 1.9484708982737828, 4.372241314917588e-05, 0.9979078845964814)

The fit can be done also using a custom optimizer:

>>> from scipy.optimize import minimize
>>> def custom_optimizer(func, x0, args=(), disp=0):
...     res = minimize(func, x0, args, method="slsqp", options={"disp": disp})
...     if res.success:
...         return res.x
...     raise RuntimeError('optimization routine failed')
>>> a1, b1, loc1, scale1 = beta.fit(x, method="MLE", optimizer=custom_optimizer)
>>> a1, b1, loc1, scale1
(1.0198821087258905, 1.948484145914738, 4.3705304486881485e-05, 0.9979104663953395)

We can also use some prior knowledge about the dataset: let's keep
``loc`` and ``scale`` fixed:

>>> a1, b1, loc1, scale1 = beta.fit(x, floc=0, fscale=1)
>>> loc1, scale1
(0, 1)

We can also keep shape parameters fixed by using ``f``-keywords. To
keep the zero-th shape parameter ``a`` equal 1, use ``f0=1`` or,
equivalently, ``fa=1``:

>>> a1, b1, loc1, scale1 = beta.fit(x, fa=1, floc=0, fscale=1)
>>> a1
1

Not all distributions return estimates for the shape parameters.
``norm`` for example just returns estimates for location and scale:

>>> from scipy.stats import norm
>>> x = norm.rvs(a, b, size=1000, random_state=123)
>>> loc1, scale1 = norm.fit(x)
>>> loc1, scale1
(0.92087172783841631, 2.0015750750324668)
rø