2017-04-11 65 views
4

私はVoigtプロファイルを実験データにフィッティングするためのコードを書くのに苦労しています。今のところ合理的なフィット関数を得ることができますが、私は約1000回の自動化を行う必要があります。つまり、2回目のフィットがすべてオフになると、それはチャンスではありません。python lmfit:voigtフィッティング - out.best_fitとout.best_valuesの差

私はフィット関数をプロットするために2つの異なる方法を試しましたが、驚くべきことに私は異なる結果を得ました。まず、out.best_fitを使ってフィット関数をプロットし、out.best_values.paramを使って私のVoigt関数に最適なパラメータを与えて試しました。

あなたがここに見ることができるように:ここでは

enter image description here

は私のコードの関連する部分である:

import matplotlib.pyplot as plt 
import numpy as np 
from lmfit import Model 



absPlot = np.genfromtxt('absPlot.txt') 


def Gaussian(x, *p): 
    A, x0, sigma, zerolevel = p 
    return A*(np.sqrt(2*np.pi*sigma))*np.exp(-(x-x0)**2/(2*sigma**2))+zerolevel 

def Lorentzian(x, *p): 
    A, x0, gamma, zerolevel = p 
    return A/np.pi*gamma/((x-x0)**2+gamma**2)+zerolevel 

def VoigtConv(x, A, x0, sigma, gamma): 
    G = np.fft.fft(Gaussian(x,1,x0,sigma,0)) 
    L = np.fft.fft(Lorentzian(x,1,x0,gamma,0)) 
    V = A*np.real(np.fft.fftshift(np.fft.ifft(G*L))) 
    return V 



vModel = Model(VoigtConv) 
vModel.set_param_hint('A', value = 111, min=100, max=2*10**4) 
vModel.set_param_hint('x0', value = 1.22*10**9, min=1.20*10**9, max=1.23*10**9) 
vModel.set_param_hint('sigma', value = 4*10**8, min=3.2*10**8, max=4.4*10**8) 
vModel.set_param_hint('gamma', value = 0, min=10**4, max=1.02*10**8) 

result = vModel.fit(absPlot[50:170,1], x=absPlot[50:170,0], method='nelder') 

#calculate y-values of fit function using out.best_values.get 
yValTest = VoigtConv(absPlot[:,0], result.best_values.get('A'), result.best_values.get('x0'), result.best_values.get('sigma'), result.best_values.get('gamma'),) 

plt.figure(1) 
plt.plot(absPlot[50:170, 0],absPlot[50:170, 1], '.k', label='Measured Data') 
plt.plot(absPlot[50:170, 0], result.best_fit, linewidth=1.2, label='Voigt Fit') 
plt.plot(absPlot[50:170, 0], yValTest[50:170], linewidth=1.2, label='Voigt Fit Calc') 
legend = plt.legend(loc='upper right', fontsize=14) 
plt.show() 

編集:今、私はいくつかの実験データを追加しました。私はあなただけの結果が同一であることを確認するために2つの行を変更する必要があり、それだけで作図問題だった私の質問にコードとして

-1.681603087629120350e+09 4.637046454387655351e-04 
-1.663354269852524996e+09 -2.908692869563400066e-03 
-1.644352957839398623e+09 7.695945924723050129e-04 
-1.624863118925523758e+09 2.147814059002654818e-03 
-1.604929495271069527e+09 -8.268747736859962489e-04 
-1.584586383181550264e+09 6.721636492005151047e-04 
-1.563867753898397207e+09 -5.488994344381660419e-04 
-1.542800945031747103e+09 2.613583859942069428e-03 
-1.521409892173170567e+09 -6.122206645226549102e-04 
-1.499717791180877686e+09 -1.831306303937777933e-03 
-1.477742650837929010e+09 4.172742657008059186e-03 
-1.455503987825480700e+09 -1.077350742753465209e-03 
-1.433014528459659338e+09 -9.746246905539714253e-04 
-1.410290534347140551e+09 3.358293207787981633e-04 
-1.387343298036361217e+09 -1.721480160424810864e-03 
-1.364185915142794132e+09 7.452103826312235285e-04 
-1.340828238149460554e+09 -5.484692294885439414e-04 
-1.317279963618578434e+09 6.280543767054739580e-04 
-1.293549994341342211e+09 2.908495463458083633e-03 
-1.269645749794533491e+09 -3.553714511279682030e-04 
-1.245576345970142126e+09 6.873643906098197585e-04 
-1.221348013080484629e+09 5.367444924189098375e-03 
-1.196966453365091801e+09 4.248147950158315603e-03 
-1.172440081617766380e+09 1.928944595089532700e-03 
-1.147771342720138073e+09 1.554006057599522294e-03 
-1.122967877091434479e+09 1.351309988748720188e-04 
-1.098032985877409935e+09 -1.177569651526613964e-03 
-1.072972883079598546e+09 4.618934106922681208e-04 
-1.047791077107019663e+09 1.458312190049944714e-04 
-1.022491622713490129e+09 -9.587366131689650014e-04 
-9.970775112038934231e+08 7.305805853652050808e-04 
-9.715539874860664606e+08 5.432066323042895103e-04 
-9.459236260065821409e+08 5.100532788727692819e-04 
-9.201888013259810209e+08 1.621000633703139487e-03 
-8.943550583160078526e+08 -2.798335907886584373e-03 
-8.684219274810085297e+08 -2.958772471271814246e-03 
-8.423946560744543076e+08 6.953094501079191102e-04 
-8.162749995434625149e+08 9.466086682280482403e-04 
-7.900645735455355644e+08 -6.637143664389960444e-04 
-7.637682778761259317e+08 2.229059812149901398e-03 
-7.373840841978719234e+08 3.086449363941101438e-03 
-7.109192549794913530e+08 -1.812785475191215070e-03 
-6.843715421246361732e+08 3.739656990641191097e-03 
-6.577454675922601223e+08 1.479560706362062364e-03 
-6.310420323323129416e+08 -6.678387318561813685e-04 
-6.042621472936900854e+08 5.032547408749150041e-03 
-5.774092528726952076e+08 3.798504735577255755e-04 
-5.504832428773229122e+08 1.171747608701265283e-03 
-5.234874298304962516e+08 -2.410484407726393944e-03 
-4.964224469999290705e+08 5.062760214963339146e-03 
-4.692897397331506610e+08 4.872895908473549344e-03 
-4.420906973762660623e+08 9.920889295731436311e-04 
-4.148257717793290019e+08 -1.789046853643734992e-03 
-3.874989028483497500e+08 1.594178042153536435e-03 
-3.601077871689334512e+08 6.375396553894154689e-03 
-3.326562783036654592e+08 1.660081886774406393e-03 
-3.051446407667568922e+08 7.034939293295599215e-03 
-2.775730907297806144e+08 3.090402499282546202e-03 
-2.499444871267271638e+08 1.745101889213399515e-03 
-2.222580747788251340e+08 4.768689476823266202e-03 
-1.945175484175896347e+08 1.447527063274924681e-03 
-1.667211811270050108e+08 1.237931188785345920e-03 
-1.388699032644901276e+08 1.028169158173280031e-03 
-1.109655202580765486e+08 1.888992310516210546e-03 
-8.300981565188933909e+07 1.061218951305126224e-02 
-5.500273516961731762e+07 1.054481962650023232e-02 
-2.694601176637900621e+07 1.159739026195029735e-02 
1.162286886406528763e+06 7.460733942283327529e-03 
2.931957016280659288e+07 9.019938720957243472e-03 
5.752325775145839900e+07 9.584925115186739702e-03 
8.577627684953397512e+07 1.136006938571703685e-02 
1.140742578420177251e+08 1.114428945159785973e-02 
1.424201814249596000e+08 9.505982102907987313e-03 
1.708115470576421320e+08 1.344830138010994457e-02 
1.992467831876419783e+08 1.187327205383587395e-02 
2.277271012588270903e+08 2.076319324572744457e-02 
2.562500412377134264e+08 2.449195242208325463e-02 
2.848159275744779706e+08 2.735082151598942565e-02 
3.134232506803109050e+08 2.774201389940251020e-02 
3.420742272177910209e+08 3.827515135158365833e-02 
3.707655220897736549e+08 4.529881503233937345e-02 
3.994984570572760105e+08 4.527952171227181410e-02 
4.282715789054992795e+08 5.607927491155496186e-02 
4.570825080519530177e+08 7.082672325946916259e-02 
4.859344808534016013e+08 7.329837060720800768e-02 
5.148232715544611812e+08 9.395665439824339715e-02 
5.437521506074851751e+08 1.097432097495898429e-01 
5.727187839551054239e+08 1.278122001625859316e-01 
6.017208387353657484e+08 1.457860184384442703e-01 
6.307616330499908924e+08 1.699138073449377728e-01 
6.598379144836143255e+08 1.930198062127272129e-01 
6.889511463229213953e+08 2.220503753486034737e-01 
7.180980850745809078e+08 2.526183222255706240e-01 
7.472821066223816872e+08 2.732252963526957124e-01 
7.764999822610137463e+08 3.132072908233128894e-01 
8.057532209395858049e+08 3.599301219050624057e-01 
8.350395505187016726e+08 3.723550635944700149e-01 
8.643595502498313189e+08 4.008270680437454603e-01 
8.937128598866268396e+08 4.628457195953887826e-01 
9.230991268193773031e+08 5.020080533129509526e-01 
9.525189598886506557e+08 5.071178382056423795e-01 
9.819710689076117277e+08 5.650011397923283551e-01 
1.011453211112567186e+09 5.762678851770945965e-01 
1.040968883829161644e+09 6.018926944100451149e-01 
1.070514901516727448e+09 6.566622250024226615e-01 
1.100091908706279755e+09 6.339041992820079185e-01 
1.129699600179332495e+09 7.038960728812075907e-01 
1.159338637641995907e+09 6.882550288044064768e-01 
1.189006806825993299e+09 6.752937525890210235e-01 
1.218704779950688362e+09 6.987419909943954899e-01 
1.248432274795983315e+09 6.907257752441064991e-01 
1.278189014494000912e+09 6.994850101135580145e-01 
1.307975693067772388e+09 6.836952966775879936e-01 
1.337790113470141411e+09 6.568788751216336763e-01 
1.367632978777440071e+09 6.368623338057169958e-01 
1.7173204e+09 6.049158734230339896e-01 
1.427403990582278967e+09 5.811201617354637694e-01 
1.457329698505902767e+09 5.569073473398759022e-01 
1.487284791765115499e+09 5.321053430579075760e-01 
1.517267092710363150e+09 4.729809914007321314e-01 
1.547275393702419758e+09 4.370085751385013872e-01 
1.577310435200158834e+09 4.232518789847144469e-01 
1.607372965113873959e+09 3.794712628494046891e-01 
1.637460812329833269e+09 3.353771063047700784e-01 
1.667573755299253464e+09 3.180599728942895554e-01 
1.697715484648065090e+09 2.796006351558285030e-01 
1.727880906450053930e+09 2.570678302712340879e-01 
1.758071765062517881e+09 2.214250934036604279e-01 
1.788287855720862627e+09 1.822474465918179909e-01 
1.818528977006859779e+09 1.689963242492904527e-01 
1.848794930768650770e+09 1.524767094429279046e-01 
1.879085522043238163e+09 1.261441533600867748e-01 
1.909401541748049259e+09 1.076852276655184126e-01 
1.939741819874165773e+09 9.167076441538968279e-02 
1.970104201912080765e+09 7.592326940626355214e-02 
2.000492440672926903e+09 6.390716920624482655e-02 
2.030903403075877905e+09 5.127758746807782597e-02 
2.061337896118272543e+09 5.026416482472045866e-02 
2.091796732556270123e+09 3.841239668649201050e-02 
2.122277766411078691e+09 3.436231750283819802e-02 
2.152782803790335655e+09 2.840162401904743408e-02 
2.183309700841841698e+09 2.428993572016554733e-02 
2.213858290976141930e+09 1.951114829537704820e-02 
2.244431383255833626e+09 1.859592448911591769e-02 
2.275026838868152142e+09 1.224343933033050086e-02 
2.305642515521794319e+09 1.188516124488762753e-02 
2.336282226802390575e+09 1.058636588047172367e-02 
2.366943836343531609e+09 1.113232016110147701e-02 
2.397625204276344299e+09 9.421615444672681861e-03 
2.428329163902065277e+09 5.280895411197075381e-03 
2.459056566942729473e+09 1.084499276586590316e-02 
2.489802289655231476e+09 1.078386352703050721e-02 
2.520570173068265438e+09 5.824088483273042946e-03 
2.551359080776900768e+09 2.315018983550010383e-03 
2.582169873054609776e+09 7.239072020225852284e-03 
2.613000416498261929e+09 6.630457173037925325e-03 
2.643850574351409435e+09 1.958339596773853771e-03 
2.674724214649637222e+09 4.983436552408344704e-03 
2.705616205165535450e+09 1.029229546263580371e-03 
2.736528416622229099e+09 7.705997351215921817e-03 
2.767460721287646770e+09 2.578113083849333659e-03 
2.798414999787809849e+09 2.746086081272369056e-03 
2.829386111613266945e+09 5.270105868119987109e-03 
2.860376942053360462e+09 1.365007411467296112e-04 
2.891390385640451908e+09 -1.043731033794107919e-03 
2.922421298921163082e+09 3.604316309236057971e-03 
2.953471572954289436e+09 5.755712938524244716e-04 
2.984540083907854557e+09 -1.467985991809585516e-04 
3.015630745689946175e+09 3.115355222446759137e-03 
3.046738408209151268e+09 1.132141324233701006e-03 
3.077864973329495907e+09 1.245146002015034333e-03 
3.109011339613405704e+09 -3.622397431994643906e-03 
3.140175378833959103e+09 7.821255833579880516e-04 
3.171357991632909775e+09 -5.220421007009016026e-03 
3.202558059600130081e+09 4.608851918083738878e-04 
3.233778511142369270e+09 -4.659729015617194645e-03 
3.265017220194062233e+09 -7.142070605820295334e-04 
3.296273070348468304e+09 -1.764112856771348501e-03 
3.327546971555518150e+09 1.964762880795126610e-03 
3.358839837300690174e+09 -1.179418031405343173e-03 
3.390148524978222370e+09 -2.105109222541607077e-03 
3.421476994503089428e+09 1.864912698552138046e-04 
3.452823119865895271e+09 -6.242222955266879016e-03 
3.484185788200661182e+09 1.285998215971130630e-04 
3.515565919819305897e+09 -3.990243330170400657e-03 
3.546963420644512177e+09 5.210052127076624136e-04 
3.578378197658686161e+09 -6.989012079038983867e-03 
3.609810158886227608e+09 -4.758719542564509609e-03 
3.641259213376206875e+09 -3.850901340059344723e-03 
3.672727311078483105e+09 3.148179067865766474e-03 
3.704210284348176003e+09 -2.667435920221450808e-03 
3.735709062965043545e+09 6.364022737269202498e-05 
3.767225601444105625e+09 1.447985623477719263e-03 
3.798760837450304031e+09 -4.522403590701212263e-03 
3.830309575790696621e+09 -2.575804324253561399e-04 
3.861874796974990845e+09 -3.712782930144581991e-03 
3.893460512536677837e+09 -3.738769907892062604e-03 
3.925058452424866199e+09 -1.397493963458843825e-03 
3.956673652574840069e+09 -1.804394360604909394e-03 
3.988305008549614906e+09 -1.262065037270537933e-03 
4.019953466279759884e+09 4.581330561478126791e-03 
4.051615869289734840e+09 -2.210331637481825415e-03 
4.083295216776277065e+09 -3.297599342719983052e-03 
4.114992460210449696e+09 1.670171676951545464e-04 
4.146704442429268837e+09 -4.767725489321614912e-03 
4.178430057674149990e+09 -2.079767584120280084e-03 
4.210171285574518204e+09 -2.953665797634475679e-03 
4.241932169318356514e+09 9.852334937458929846e-04 
4.273704402316565514e+09 -1.863814866179199533e-04 
4.305492027227878571e+09 -1.995553593671882224e-03 
4.337298064545919418e+09 -4.032630887044518289e-03 
4.369117290490421295e+09 3.192997791401750059e-03 
4.400950663249242783e+09 -4.609954747607243782e-03 
4.432801208367661476e+09 -6.043342046874808230e-04 
4.464667826607993126e+09 -9.782805973301950640e-04 
4.496546317504646301e+09 -2.979962925485028519e-03 
4.528440743885611534e+09 -6.625228674059118784e-04 
4.560350005576179504e+09 -2.660491280681023633e-03 
4.592275070624813080e+09 -1.162824608277974125e-04 
4.624212768708718300e+09 -1.468441922582553122e-03 
4.656166138638625145e+09 -4.418446727839028115e-03 
4.688136153178585052e+09 -8.467118791638769443e-05 
4.720118604362720490e+09 2.264458645799919626e-03 
4.752114462863739014e+09 -4.646356635505886879e-04 
4.784126776974655151e+09 5.455519464959710310e-04 
4.816154448997751236e+09 -1.487020291985546557e-03 
4.848192225956024170e+09 3.737703214782029607e-03 
4.880246274756915092e+09 3.274430116837369256e-03 
4.912316537408501625e+09 -4.280346194170173840e-03 
4.944398796976040840e+09 -1.503316968807124135e-03 
4.976494031398631096e+09 6.984472814897352104e-04 
+0

あなたはまた、データを提供することができ、それ以外の場合はあなたの問題を再現するのは難しいです。 – Cleb

+1

私はちょうどいくつかのデータを追加しました。私はコードブロックとして追加するよりも良いアイデアはありませんでした。 – oigen

+0

私は 'result.best_fit'と' yValTest'の結果を見て、私が知る限り同じものでした。どちらの場合も、私は '[0.00340746、0.00345166、0.00351154 ... 0.0033669、0.00336503、0.00337859] 'を受け取っています(最初と最後の3つの値のみを表示します)。それではまったく違うのは何ですか? :) – Cleb

答えて

4

それをコピーして貼り付けるよりも、より良いアイデアを持っていませんでした:

yValTest = VoigtConv(absPlot[:,0], result.best_values.get('A'), result.best_values.get('x0'), result.best_values.get('sigma'), result.best_values.get('gamma'),) 

またその後

yValTest = VoigtConv(absPlot[50:170, 0], result.best_values.get('A'), result.best_values.get('x0'), result.best_values.get('sigma'), result.best_values.get('gamma'),) 

(変更されたスライスに注意)及び

plt.plot(absPlot[50:170, 0], yValTest[50:170], linewidth=1.2, label='Voigt Fit Calc') 

plt.plot(absPlot[50:170, 0], yValTest, linewidth=1.2, label='Voigt Fit Calc') 

(変更されたスライスに注意してください)に問題はあなたがお互いに異なる範囲を比較したということでした。これを固定した後、予想通り、プロットは(青と緑が区別できない。私は、下記の別のテストを追加しました)になります。

enter image description here

もう一つの簡単なテスト:つまり

np.array_equal(result.best_fit, yValTest) 

戻りTrue、 2つの配列は同じです。

はここで全体のコードです(私はパンダを使用してデータを読み込むには、再度この行を変更することができます):

import matplotlib.pyplot as plt 
import numpy as np 
from lmfit import Model 
import pandas as pd 


absPlot = pd.read_csv('lmfit_voigt.csv', index_col=None).values 


def Gaussian(x, *p): 
    A, x0, sigma, zerolevel = p 
    return A*(np.sqrt(2*np.pi*sigma))*np.exp(-(x-x0)**2/(2*sigma**2))+zerolevel 

def Lorentzian(x, *p): 
    A, x0, gamma, zerolevel = p 
    return A/np.pi*gamma/((x-x0)**2+gamma**2)+zerolevel 

def VoigtConv(x, A, x0, sigma, gamma): 
    G = np.fft.fft(Gaussian(x,1,x0,sigma,0)) 
    L = np.fft.fft(Lorentzian(x,1,x0,gamma,0)) 
    V = A*np.real(np.fft.fftshift(np.fft.ifft(G*L))) 
    return V 


vModel = Model(VoigtConv) 
vModel.set_param_hint('A', value = 111, min=100, max=2*10**4) 
vModel.set_param_hint('x0', value = 1.22*10**9, min=1.20*10**9, max=1.23*10**9) 
vModel.set_param_hint('sigma', value = 4*10**8, min=3.2*10**8, max=4.4*10**8) 
vModel.set_param_hint('gamma', value = 0, min=10**4, max=1.02*10**8) 

result = vModel.fit(absPlot[50:170,1], x=absPlot[50:170,0], method='nelder') 

#calculate y-values of fit function using out.best_values.get 
yValTest = VoigtConv(absPlot[50:170, 0], result.best_values.get('A'), result.best_values.get('x0'), result.best_values.get('sigma'), result.best_values.get('gamma'),) 

plt.figure(1) 
plt.plot(absPlot[50:170, 0], absPlot[50:170, 1], '.k', label='Measured Data') 
plt.plot(absPlot[50:170, 0], result.best_fit, linewidth=1.2, label='Voigt Fit') 
plt.plot(absPlot[50:170, 0], yValTest, linewidth=1.2, label='Voigt Fit Calc') 
legend = plt.legend(loc='upper right', fontsize=14) 
plt.show() 
+1

Cleb:とても素敵な答えに感謝します! –

+1

@Clebこの詳細な回答をありがとう!私は同様のプロットを取得し、私はもはや混乱していないことを知っている。 – oigen

関連する問題