scipy powerlaw fit에서 합리적인 값을 얻으려고합니다.

Question

전 전력 법 의존도를 계산하기 위해 실행중인 시뮬레이션 코드의 일부 데이터를 맞추려고합니다. 선형 적합을 플롯 할 때 데이터가 잘 맞지 않습니다.

#!/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

반면에 여기

내가 데이터에 맞게 사용하고 파이썬 스크립트입니다 줄거리가 최소한으로 잘 맞는 것으로 기대되는만큼 잘 맞는 것처럼 보이고, 출력의 형태가 저를 귀찮게합니다. 나는 상수가 당신이 0이 될 것으로 기대하는 곳에 가깝기를 바랐다 (약 30). 그리고 나는 10^-5보다 더 큰 부분의 힘 의존성을 발견 할 것을 기대하고 있었다.

나는 데이터를 재조정하고 매개 변수를 가지고 놀아서 행운이없는 optimize.leastsq를 시도했다. 가능한 목표를 달성하려고합니까 아니면 내 데이터가 허용하지 않습니까? 계산이 비싸기 때문에 더 많은 데이터 포인트를 얻는 것이 중요하지 않습니다.

감사합니다.

Answer 1

숫자가 모두 너무 작지 않도록 xdata의 크기를 조정하면 도움이됩니다. 새 변수 xprime = 1000*x에서 작업 할 수 있습니다. 다음 xprime 대 y을 맞습니다.

최소 제곱 그래서 같은 그런

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

또한 원하는 결과에 가까운 뭔가 초기 추측을 변경하는 데 도움이 y = p[0] + p[1] * (x ** p[2])

,하자

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

피팅 매개 변수 q를 찾을 수 [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

Answer 2

먼저, 대수를 취하면 훨씬를 제공이 차식에 맞게 leastsquare을 사용하는 것이 훨씬 더 낫다 더 잘 어울린다. scipy cookbook에는 훌륭한 예가 있습니다. 아래에서 귀하의 코드에 맞도록 수정했습니다.

같은 가장 맞는

은 우리가 그래프 (및 데이터)에서 볼 수 있듯이 진폭 = 0.8955, 및 인덱스 = -0.40943265484

, 그 경우 전원 법에 맞게 우리는 진폭 값을 기대하지 않을 것이다 30 근처에 있습니다. 멱 법칙 방정식 f(x) == Amp * x ** index 에서처럼 음수 지수는 f(1) == Amp 및 f(0) == infinity입니다. 로그인 할 때 (y_i)를 직선에 맞게 :

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() 

Answer 3

지수 적합을 얻기 위해 선형 최소 제곱을 사용하는 표준 방법은 무엇 fraxel suggests in his/her answer하는 것입니다.

그러나이 방법은 수치적인 단점, 특히 감도 (데이터의 작은 변화로 인해 견적에서 큰 변화를 가져옴)를 가지고 있습니다. 선호되는 대안은 비선형 최소 제곱 접근법을 사용하는 것입니다. 이는 덜 민감합니다. 그러나 중요하지 않은 목적으로 선형 LS 방법에 만족한다면 그냥 사용하십시오.