Пытаясь получить разумные значения от scipy powerlaw

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

Здесь python script, который я использую, чтобы соответствовать данным:

#!/usr/bin/env python
from scipy import optimize
import numpy

xdata=[ 0.00010851,  0.00021701,  0.00043403,  0.00086806,  0.00173611, 0.00347222]
ydata=[ 29.56241016,  29.82245508,  25.33930469,  19.97075977,  12.61276074, 7.12695312]

fitfunc = lambda p, x: p[0] + p[1] * x ** (p[2])
errfunc = lambda p, x, y: (y - fitfunc(p, x))

out,success = optimize.leastsq(errfunc, [1,-1,-0.5],args=(xdata, ydata),maxfev=3000)

print "%g + %g*x^%g"%(out[0],out[1],out[2])

вывод, который я получаю: -71205.3 + 71174.5 * x ^ -9.79038e-05

В то время как на сюжете пригонка выглядит так же хорошо, как вы ожидали бы от наименьших квадратов, форма вывода беспокоит меня. Я надеялся, что константа будет близка к тому, где вы ожидаете, что ноль будет (около 30). И я ожидал найти зависимость мощности большей доли, чем 10 ^ -5.

Я пробовал перемасштабировать свои данные и играть с параметрами optimize.leastsq без везения. Является ли то, что я пытаюсь сделать возможным, или мои данные просто не позволяют? Расчет дорог, поэтому получение большего количества точек данных является нетривиальным.

Спасибо!

Ответ 1

Это помогает перемасштабировать xdata, поэтому числа не так малы. Вы можете работать с новой переменной xprime = 1000*x. Затем установите xprime в сравнении с y.

Наименьшие квадраты найдут параметры q подгонка

y = q[0] + q[1] * (xprime ** q[2]) 
  = q[0] + q[1] * ((1000*x) ** q[2])

Итак, пусть

p[0] = q[0]
p[1] = q[1] * (1000**q[2])
p[2] = q[2]

Тогда y = p[0] + p[1] * (x ** p[2])

Это также помогает изменить первоначальное предположение на нечто более близкое к вашему желаемому результату, например [max(ydata), -1, -0.5].

from scipy import optimize
import numpy as np

def fitfunc(p, x):
    return p[0] + p[1] * (x ** p[2])
def errfunc(p, x, y):
    return y - fitfunc(p, x)

xdata=np.array([ 0.00010851,  0.00021701,  0.00043403,  0.00086806,
                 0.00173611, 0.00347222])
ydata=np.array([ 29.56241016,  29.82245508,  25.33930469,  19.97075977,
                 12.61276074, 7.12695312])

N = 5000
xprime = xdata * N

qout,success = optimize.leastsq(errfunc, [max(ydata),-1,-0.5],
                               args=(xprime, ydata),maxfev=3000)

out = qout[:]
out[0] = qout[0]
out[1] = qout[1] * (N**qout[2])
out[2] = qout[2]
print "%g + %g*x^%g"%(out[0],out[1],out[2])

дает

40.1253 + -282.949 * x ^ 0.375555

Ответ 2

Лучше сначала взять логарифм, а затем использовать leastsquare, чтобы соответствовать этому линейному уравнению, что даст вам гораздо лучшую форму. В scipy cookbook есть отличный пример, который я адаптировал ниже, чтобы соответствовать вашему коду.

Наиболее подходящими являются такие: амплитуда = 0,8955, а индекс = -0.40943265484

Как видно из графика (и ваших данных), если его степенной закон подходит, мы не ожидаем, что значение амплитуды будет близким к 30. Как и в уравнении степенного закона f(x) == Amp * x ** index, поэтому с отрицательным индексом: f(1) == Amp и f(0) == infinity.

from pylab import *
from scipy import *
from scipy import optimize

xdata=[ 0.00010851,  0.00021701,  0.00043403,  0.00086806,  0.00173611, 0.00347222]
ydata=[ 29.56241016,  29.82245508,  25.33930469,  19.97075977,  12.61276074, 7.12695312]

logx = log10(xdata)
logy = log10(ydata)

# define our (line) fitting function
fitfunc = lambda p, x: p[0] + p[1] * x
errfunc = lambda p, x, y: (y - fitfunc(p, x))

pinit = [1.0, -1.0]
out = optimize.leastsq(errfunc, pinit,
                       args=(logx, logy), full_output=1)

pfinal = out[0]
covar = out[1]

index = pfinal[1]
amp = 10.0**pfinal[0]

print 'amp:',amp, 'index', index

powerlaw = lambda x, amp, index: amp * (x**index)
##########
# Plotting data
##########
clf()
subplot(2, 1, 1)
plot(xdata, powerlaw(xdata, amp, index))     # Fit
plot(xdata, ydata)#, yerr=yerr, fmt='k.')  # Data
text(0.0020, 30, 'Ampli = %5.2f' % amp)
text(0.0020, 25, 'Index = %5.2f' % index)
xlabel('X')
ylabel('Y')

subplot(2, 1, 2)
loglog(xdata, powerlaw(xdata, amp, index))
plot(xdata, ydata)#, yerr=yerr, fmt='k.')  # Data
xlabel('X (log scale)')
ylabel('Y (log scale)')

savefig('power_law_fit.png')
show()

Ответ 3

Стандартный способ использования линейных наименьших квадратов для получения экспоненциального соответствия заключается в том, чтобы сделать то, что fraxel предлагает в своем ответе: поместить прямую линию в журнал (y_i).

Однако этот метод имеет известные численные недостатки, особенно чувствительность (небольшое изменение в данных дает большое изменение в оценке). Предпочтительным вариантом является использование нелинейного метода наименьших квадратов - он менее чувствителен. Но если вы удовлетворены линейным методом LS для некритических целей, просто используйте это.