Как настроить функцию из `scipy.stats` для работы с` scipy.optimize.curve_fit`? - PullRequest
0 голосов
/ 04 мая 2018

Я создаю единый вектор вероятностей, добавляю вес к региону, снова превращаюсь в вероятности. Я хочу сгладить кривую, подгоняя распределение beta к кривой. Я не могу использовать stats.beta.fit, потому что у меня нет ничьих из этих вероятностей, только каркас для кривой.

Как мне настроить функцию из scipy.stats для работы с scipy.optimize.curve_fit? Я не хочу ограничивать ее только beta дистрибутивами, если это возможно. Есть ли общий способ преобразовать их в форму, которую можно минимизировать для оптимизации набора параметров?

По сути, я ищу набор параметров, который наилучшим образом соответствует кривой для конкретного распределения (в данном случае beta).

Есть идеи?

# Uniform probabilities
n = 10
x = np.linspace(0,1,n)
probs_num = np.ones(n)/n # y-values

# Update the mass around an area
idx_update = 2
probs_num[idx_update] += 0.382

# Convert fo probs
probs_num = probs_num/probs_num.sum()

# Plot
with plt.style.context("ggplot"):
    fig, ax = plt.subplots()
    ax.plot(x,probs_num)
    ax.set_ylabel("Density")
    ax.set_xlabel("$x$")


probs_num

enter image description here

from scipy.optimize import curve_fit
from scipy import stats
popt, pcov = curve_fit(stats.beta, x, probs_num) # NOTE: I know this is the wrong way to use it but I'm leaving it as a placeholder

В ответ на ответ lm ниже:

import lmfit
def beta_fcn(x, alpha, beta, loc):
    return stats.beta.pdf(x, alpha, beta, loc)
bmodel = lmfit.Model(beta_fcn)
params = bmodel.make_params(alpha=1, beta=1., loc=0.5)
result = bmodel.fit(probs_num, params, x=x)

# ---------------------------------------------------------------------------
# ValueError                                Traceback (most recent call last)
# <ipython-input-34-61ec32095934> in <module>()
#       4 bmodel = lmfit.Model(beta_fcn)
#       5 params = bmodel.make_params(alpha=1, beta=1., loc=0.5)
# ----> 6 result = bmodel.fit(probs_num, params, x=x)

# ~/anaconda/envs/python3/lib/python3.6/site-packages/lmfit/model.py in fit(self, data, params, weights, method, iter_cb, scale_covar, verbose, fit_kws, nan_policy, **kwargs)
#     871                              scale_covar=scale_covar, fcn_kws=kwargs,
#     872                              nan_policy=self.nan_policy, **fit_kws)
# --> 873         output.fit(data=data, weights=weights)
#     874         output.components = self.components
#     875         return output

# ~/anaconda/envs/python3/lib/python3.6/site-packages/lmfit/model.py in fit(self, data, params, weights, method, nan_policy, **kwargs)
#    1215         self.userkws.update(kwargs)
#    1216         self.init_fit = self.model.eval(params=self.params, **self.userkws)
# -> 1217         _ret = self.minimize(method=self.method)
#    1218 
#    1219         for attr in dir(_ret):

# ~/anaconda/envs/python3/lib/python3.6/site-packages/lmfit/minimizer.py in minimize(self, method, params, **kws)
#    1636                         val.lower().startswith(user_method)):
#    1637                     kwargs['method'] = val
# -> 1638         return function(**kwargs)
#    1639 
#    1640 

# ~/anaconda/envs/python3/lib/python3.6/site-packages/lmfit/minimizer.py in leastsq(self, params, **kws)
#    1288         np.seterr(all='ignore')
#    1289 
# -> 1290         lsout = scipy_leastsq(self.__residual, variables, **lskws)
#    1291         _best, _cov, infodict, errmsg, ier = lsout
#    1292         result.aborted = self._abort

# ~/anaconda/envs/python3/lib/python3.6/site-packages/scipy/optimize/minpack.py in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
#     385             maxfev = 200*(n + 1)
#     386         retval = _minpack._lmdif(func, x0, args, full_output, ftol, xtol,
# --> 387                                  gtol, maxfev, epsfcn, factor, diag)
#     388     else:
#     389         if col_deriv:

# ~/anaconda/envs/python3/lib/python3.6/site-packages/lmfit/minimizer.py in __residual(self, fvars, apply_bounds_transformation)
#     489         if not self._abort:
#     490             return _nan_policy(np.asarray(out).ravel(),
# --> 491                                nan_policy=self.nan_policy)
#     492 
#     493     def __jacobian(self, fvars):

# ~/anaconda/envs/python3/lib/python3.6/site-packages/lmfit/minimizer.py in _nan_policy(arr, nan_policy, handle_inf)
#    1846 
#    1847         if contains_nan:
# -> 1848             raise ValueError("The input contains nan values")
#    1849     return arr
#    1850 

# ValueError: The input contains nan values

1 Ответ

0 голосов
/ 05 мая 2018

Возможно, вы захотите использовать beta.pdf, например:

import numpy as np
from scipy import stats
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

n = 21
x = np.linspace(0, 1, n)
y = stats.beta.pdf(x, 2.5, 15, 0.3) + np.random.normal(scale=0.1, size=n)

p0 = [2.0, 10.0, 0.5]
popt, pcov = curve_fit(stats.beta.pdf, x, y, p0)
print(popt)

Вы также можете найти lmfit (https://lmfit.github.io/lmfit-py/) полезным. Для этого вам нужно будет обернуть stats.beta.pdf, но это дает вам гораздо большую гибкость в установлении границ параметров или их исправлении:

from lmfit import Model
def beta_fcn(x, alpha, beta, loc):
    return stats.beta.pdf(x, alpha, beta, loc)

bmodel = Model(beta_fcn)

params = bmodel.make_params(alpha=2, beta=10., loc=0.5)
params['alpha'].min = 1.0

result = bmodel.fit(y, params, x=x)

print(result.fit_report())

with plt.style.context("ggplot"):
    fig, ax = plt.subplots()
    ax.plot(x, y, 'o')
    ax.plot(x, result.best_fit)
    ax.set_ylabel("Density")
    ax.set_xlabel("$x$")
plt.show()

Это распечатает

[[Model]]
    Model(beta_fcn)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 75
    # data points      = 21
    # variables        = 3
    chi-square         = 0.06244742
    reduced chi-square = 0.00346930
    Akaike info crit   = -116.177008
    Bayesian info crit = -113.043441
[[Variables]]
    alpha:  2.50445306 +/- 0.06843391 (2.73%) (init = 2)
    beta:   14.9105872 +/- 0.31198655 (2.09%) (init = 10)
    loc:    0.29906057 +/- 0.00190835 (0.64%) (init = 0.5)
[[Correlations]] (unreported correlations are < 0.100)
    C(alpha, beta) =  0.894
    C(alpha, loc)  = -0.834
    C(beta, loc)   = -0.564

(и, чтобы быть ясным, наилучшие значения соответствия и неопределенности будут такими же, как в примере curve_fit выше) и график

enter image description here

...