scipyのダウンロードカーブフィットだからべき乗則

Question

を失敗し、私は次のような種類のべき乗則でデータのセットに合うようにしようとしている：

def f(x,N,a): # Power law fit 
    if a >0: 
     return N*x**(-a) 
    else: 
     return 10.**300 

par,cov = scipy.optimize.curve_fit(f,data,time,array([10**(-7),1.2])) 

他の条件はちょうど正であることを強制することです。 scipy.optimize.curve_fitを使用すると、データとの交差点がまったくなく、それぞれNとaに対して1.2e + 04と1.9e0-7の値を返すan awful fit (green line)が得られます。私が手動で入力したフィットから、Nとaの値はそれぞれ1e-07と1.2になるはずですが、初期パラメータとしてcurve_fitに入れても結果は変わりません。 aが正であるという条件を削除すると、負の値が選択されるため、フィットが悪くなります。これは、間違った符号の傾きでフィットすることにつながります。

私は信じられないほど信頼できるものの、このルーチンから外れないようにする方法を理解することはできませんが、他のPythonカーブフィッティングルーチンを見つけることはできません。自分の最小自乗アルゴリズムを書く必要があるのですか、ここで間違っていることがありますか？オリジナルのポストで

Answer 1

UPDATE

、私はあなたのパラメータに境界を割り当てることができますlmfitを使用するソリューションを示しました。バージョン0.17以降、scipyはパラメータに直接境界を割り当てることもできます（documentation参照）。 EDITの後ろにあるこのソリューションは、scipyのcurve_fitとパラメータ境界の使用方法に関する最小限の例として役立ちます。

オリジナルポスト

@Warren Weckesserによって示唆されるように、あなたはあなたのパラメータに境界を割り当てることができますし、この「醜い」IF-句を回避する、このタスクを成し遂げるためにlmfitを使用することができます。

あなたが任意のデータを提供していないので、私はここに示されているいくつかの作成：彼らは私が彼らに合わせ、法律にf(x) = 10.5 * x ** (-0.08)

に従ってください

を変換することによって - - @のroadrunner66により示唆されるように線形関数でのべき乗則：

y = N * x ** a 
ln(y) = ln(N * x ** a) 
ln(y) = a * ln(x) + ln(N) 

だから私は、最初のoriにnp.logを使用ginal dataと入力してからフィットを行います。私は今lmfitを使用すると、私は次のような出力が得られます。

[[Variables]] 
    lN: 2.35450302 +/- 0.019531 (0.83%) (init= 1.704748) 
    a: -0.08035342 +/- 0.005158 (6.42%) (init=-0.5) 

だからaは元の値にかなり近く、np.exp(2.35450302)も元の値に非常に近い10.53を与えます。

プロットは次のようになります。あなたが見ることができるようにフィットは非常によくデータを示しています。ここでは

は、インラインコメントのカップルと全体のコードです：

import numpy as np 
import matplotlib.pyplot as plt 
from lmfit import minimize, Parameters, Parameter, report_fit 

# generate some data with noise 
xData = np.linspace(0.01, 100., 50.) 
aOrg = 0.08 
Norg = 10.5 
yData = Norg * xData ** (-aOrg) + np.random.normal(0, 0.5, len(xData)) 
plt.plot(xData, yData, 'bo') 
plt.show() 

# transform data so that we can use a linear fit 
lx = np.log(xData) 
ly = np.log(yData) 
plt.plot(lx, ly, 'bo') 
plt.show() 

def decay(params, x, data): 

    lN = params['lN'].value 
    a = params['a'].value 

    # our linear model 
    model = a * x + lN 
    return model - data # that's what you want to minimize 

# create a set of Parameters 
params = Parameters() 
params.add('lN', value=np.log(5.5), min=0.01, max=100) # value is the initial value 
params.add('a', value=-0.5, min=-1, max=-0.001) # min, max define parameter bounds 

# do fit, here with leastsq model 
result = minimize(decay, params, args=(lx, ly)) 

# write error report 
report_fit(params) 

# plot data 
xnew = np.linspace(0., 100., 5000.) 
# plot the data 
plt.plot(xData, yData, 'bo') 
plt.plot(xnew, np.exp(result.values['lN']) * xnew ** (result.values['a']), 'r') 
plt.show() 

EDIT

は、あなたと仮定すると、 scipy 0.17がインストールされている場合は、curve_fitを使用して次の操作を行うこともできます。私は対数データ（下のプロットの黒い線）と同様に、べき乗法（下のプロットの赤線）の元の定義のためにそれを示します。データは上記と同じ方法で生成されます。プロットとしてルックスは次のとおりです。

あなたが見ることができるように、データは非常によく説明しています。 poptとpopt_logを印刷すると、それぞれarray([ 10.47463426, 0.07914812])とarray([ 2.35158653, -0.08045776])が得られます（注：文字1の場合、最初の引数の指数をとる必要があります - np.exp(popt_log[0]) = 10.502は元のデータに近いです）。ここで

は、全体のコードです：

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

# generate some data with noise 
xData = np.linspace(0.01, 100., 50) 
aOrg = 0.08 
Norg = 10.5 
yData = Norg * xData ** (-aOrg) + np.random.normal(0, 0.5, len(xData)) 

# get logarithmic data 
lx = np.log(xData) 
ly = np.log(yData) 

def f(x, N, a): 
    return N * x ** (-a) 

def f_log(x, lN, a): 
    return a * x + lN 

# optimize using the appropriate bounds 
popt, pcov = curve_fit(f, xData, yData, bounds=(0, [30., 20.])) 
popt_log, pcov_log = curve_fit(f_log, lx, ly, bounds=([0, -10], [30., 20.])) 

xnew = np.linspace(0.01, 100., 5000) 

# plot the data 
plt.plot(xData, yData, 'bo') 
plt.plot(xnew, f(xnew, *popt), 'r') 
plt.plot(xnew, f(xnew, np.exp(popt_log[0]), -popt_log[1]), 'k') 
plt.show()