Подходящий степенной закон с использованием scipy минимизировать

Question

Я хочу сделать оценку максимального правдоподобия, используя scipy minimize для функции формы

y = a * x^b

где ошибки предполагаются нормально распределенными по прогнозам. В настоящее время истинные параметры не могут быть восстановлены вообще (см. Ниже).

Функция отрицательного логарифмического правдоподобия:

def loglik(param, data):
    mu = 0.0
    logLikelihood = 0.0
    resid = 0.0
    for i in range(data.shape[0]):
        mu = param[0] * (data[i][0] ** param[1])
        resid = float(data[i][1] - mu) / param[2]
        logLikelihood += -0.5 * resid * resid

    return -(-data.shape[0] * np.log(param[2]) + logLikelihood)

Вот код, который я использую до сих пор, но, как вы увидите, параметры после конвергенции далеки от истинных.

def generate(param, x):
    pred = [(param[0] * (x ** param[1])) for x in x]
    return np.array([sum(x) for x in zip(pred, np.random.normal(0, param[2], len(x)))])

генерировать поддельные данные

x = np.linspace(1, 50, num=100)
true_param = [15, 1.1, 5]
data = np.vstack((x, generate(true_param, x))).T

подходящая модель

from scipy.optimize import minimize
initParams = [1,1,1]
result = minimize(loglik, method='Nelder-Mead', x0=initParams,args=data)

print(result.x) # should be 15, 1.1, 5