import numpy as np
import matplotlib.pyplot as plt

#Breit-Wignerovo rozdělení
def lorentzian(x,x0,l):
    return(1/np.pi*l/(l**2+(x-x0)**2))

#Věrohodnostní funkce
def logLikelihood(par,data):
    Ndata=np.size(data)
    sum=0
    for i in range(Ndata):
        sum=sum+np.log(lorentzian(data[i],par[0],par[1]))
    return(sum)    

#Výpočet Jacobiho matice
def Jacobi(par,data):
    Ndata=np.size(data)
    sum=np.zeros(2)
    for i in range(Ndata):
        sum[0]=sum[0]+2*(data[i]-par[0])/(par[1]**2+(data[i]-par[0])**2)
        sum[1]=sum[1]+1/par[1]-2*par[1]/(par[1]**2+(data[i]-par[0])**2)
    return(sum)    

#Výpočet Hesseovi matice
def Hesse(par,data):
    Ndata=np.size(data)
    sum=np.zeros([2,2])
    for i in range(Ndata):
        sum[0,0]=sum[0,0]-2*(par[1]**2-(data[i]-par[0])**2)/(par[1]**2+(data[i]-par[0])**2)**2
        sum[1,1]=sum[1,1]-1/par[1]**2+2*(par[1]**2-(data[i]-par[0])**2)/(par[1]**2+(data[i]-par[0])**2)**2
        sum[0,1]=sum[0,1]-4*par[1]*(data[i]-par[0])/(par[1]**2+(data[i]-par[0])**2)**2
    sum[1,0]=sum[0,1]    
    return(sum)

N=100
x0true=4
ltrue=1

#simulace
x=np.empty(N)
y_point=np.zeros(N)+1/(np.pi*ltrue)
x=ltrue*np.tan((np.random.uniform(size=N)-0.5)*np.pi)+x0true
x_plot=np.linspace(x0true-10*ltrue,x0true+10*ltrue,100)
y_plot=lorentzian(x_plot,x0true,ltrue) 
fig,ax=plt.subplots()
ax.hist(x,bins=2*int(np.sqrt(N)),range=(x0true-10*ltrue,x0true+10*ltrue),density='true')
ax.scatter(x,y_point,c='black')
plt.xlim(x0true-10*ltrue,x0true+10*ltrue)
ax.plot(x_plot,y_plot,c='red')
ax.set_xlabel('x',fontsize=12)
ax.set_ylabel('četnost',fontsize=12)
plt.show()

# intial estimate
par_estimate=np.empty(2)
par_estimate[0]=4.2
par_estimate[1]=1.1
par=np.empty(2)
eps=1e-9

diff=np.sqrt(np.dot(Jacobi(par_estimate,x),Jacobi(par_estimate,x)))
fig,ax=plt.subplots(3,figsize=(8,8))
fig.tight_layout(h_pad=2)
#ax[0].set_xlabel('iterace',fontsize=12)
ax[0].set_ylabel('x0',fontsize=12)
#ax[1].set_xlabel('iterace',fontsize=12)
ax[1].set_ylabel('l',fontsize=12)
ax[2].set_xlabel('iterace',fontsize=12)
ax[2].set_ylabel('| J |',fontsize=12)
#Newton-Raphson
istep=0
while diff>eps: 
    par=par_estimate-0.1*np.matmul(np.linalg.inv(Hesse(par_estimate,x)),Jacobi(par_estimate,x))
    diff=np.sqrt(np.dot(Jacobi(par_estimate,x),Jacobi(par_estimate,x)))
    ax[0].scatter(istep,par_estimate[0],c='red')
    ax[1].scatter(istep,par_estimate[1],c='green')
    ax[2].scatter(istep,diff,c='blue')
    print(par,diff)
    par_estimate=np.copy(par)
    istep+=1
ax[0].set_title('poloha majáku x0')
ax[1].set_title('vzdálenost majáku l')
ax[2].set_title('velikost Jacobiho matice |J|')
plt.show()

[4.19732557 1.09693063] 1.6572856694689622
[4.19493532 1.09420116] 1.4908103773891956
[4.19279735 1.09177082] 1.3411314721979264
[4.19088371 1.08960435] 1.2065390426053535
[4.18916982 1.08767115] 1.0855004950585008
[4.18763399 1.08594455] 0.9766414510728602
[4.18625709 1.08440123] 0.8787288659031607
[4.18502215 1.08302076] 0.7906560717806088
[4.18391412 1.08178515] 0.7114294977073606
[4.18291964 1.08067859] 0.6401568564615184
[4.1820268  1.07968708] 0.5760366205683703
[4.18122502 1.07879827] 0.518348634346954
[4.18050483 1.07800118] 0.46644573004869705
[4.17985781 1.0772861 ] 0.4197462335304171
[4.17927641 1.07664438] 0.37772725955982694
[4.17875389 1.07606832] 0.3399187092816582
[4.17828423 1.07555107] 0.3058978929876195
[4.17786201 1.07508652] 0.2752847104596938
[4.1774824  1.07466921] 0.2477373290428886
[4.17714107 1.07429426] 0.22294830645426844
[4.17683412 1.07395731] 0.20064111130982665
[4.17655808 1.07365447] 0.18056699958155217
[4.1763098  1.07338226] 0.16250220979171778
[4.17608649 1.07313753] 0.14624544379762902
[4.17588561 1.07291749] 0.13161560359372007
[4.17570491 1.07271964] 0.11844975771970481
[4.17554235 1.07254171] 0.10660131366725976
[4.17539611 1.07238169] 0.09593837516737197
[4.17526453 1.07223777] 0.08634226545563559
[4.17514615 1.07210831] 0.07770619958471273
[4.17503964 1.07199186] 0.06993409061208025
[4.1749438 1.0718871] 0.06293947606075515
[4.17485757 1.07179287] 0.05664455245298122
[4.17477998 1.07170808] 0.05097930697070431
[4.17471015 1.0716318 ] 0.045880736418013324
[4.17464733 1.07156318] 0.0412921446653447
[4.17459079 1.07150143] 0.03716251065443673
[4.17453991 1.07144587] 0.033445919849465544
[4.17449413 1.07139587] 0.03010105274210266
[4.17445293 1.07135089] 0.02709072466726969
[4.17441585 1.07131041] 0.024381471767412167
[4.17438249 1.07127398] 0.02194317846559949
[4.17435246 1.07124121] 0.019748742276527117
[4.17432544 1.07121171] 0.017773772205128352
[4.17430112 1.07118517] 0.015996317361147303
[4.17427923 1.07116128] 0.014396622757327235
[4.17425954 1.07113979] 0.012956909564039781
[4.17424181 1.07112044] 0.011661177368273011
[4.17422586 1.07110303] 0.010495026230382737
[4.1742115  1.07108737] 0.00944549655454051
[4.17419858 1.07107327] 0.008500924987810621
[4.17418695 1.07106058] 0.007650814741970785
[4.17417649 1.07104916] 0.006885718893486094
[4.17416707 1.07103889] 0.0061971353615230705
[4.17415859 1.07102964] 0.005577412395283199
[4.17415097 1.07102131] 0.0050196635176850204
[4.1741441  1.07101382] 0.004517690979313282
[4.17413792 1.07100708] 0.004065916870405573
[4.17413236 1.07100101] 0.003659321124586858
[4.17412736 1.07099555] 0.0032933857245969846
[4.17412285 1.07099064] 0.0029640444893124234
[4.1741188  1.07098622] 0.002667637883530736
[4.17411515 1.07098224] 0.0024008723481627812
[4.17411187 1.07097865] 0.002160783698279061
[4.17410891 1.07097543] 0.0019447041822632143
[4.17410625 1.07097253] 0.0017502328356584545
[4.17410386 1.07096992] 0.001575208800096467
[4.1741017  1.07096757] 0.00141768731099002
[4.17409977 1.07096545] 0.0012759180865214545
[4.17409802 1.07096355] 0.001148325878257207
[4.17409645 1.07096184] 0.0010334929667304893
[4.17409504 1.07096029] 0.0009301434078469398
[4.17409376 1.07095891] 0.0008371288546982147
[4.17409262 1.07095766] 0.0007534157972071762
[4.17409159 1.07095653] 0.0006780740781511047
[4.17409066 1.07095552] 0.0006102665574866626
[4.17408983 1.07095461] 0.0005492398103261107
[4.17408908 1.07095379] 0.0004943157552243048
[4.1740884  1.07095305] 0.00044488411971362204
[4.17408779 1.07095239] 0.0004003956591552969
[4.17408724 1.07095179] 0.00036035605389393195
[4.17408675 1.07095126] 0.0003243204166372619
[4.17408631 1.07095077] 0.00029188834915581463
[4.17408591 1.07095034] 0.0002626994933365411
[4.17408555 1.07094994] 0.0002364295270605077
[4.17408523 1.07094959] 0.0002127865606313331
[4.17408494 1.07094927] 0.00019150789347803567
[4.17408467 1.07094899] 0.00017235709513261472
[4.17408444 1.07094873] 0.00015512137832036422
[4.17408423 1.0709485 ] 0.0001396092345709098
[4.17408403 1.07094829] 0.0001256483063276877
[4.17408386 1.0709481 ] 0.0001130834718287181
[4.17408371 1.07094794] 0.00010177512149427137
[4.17408357 1.07094778] 9.159760681779763e-05
[4.17408344 1.07094765] 8.243784405564359e-05
[4.17408333 1.07094752] 7.41940579913821e-05
[4.17408323 1.07094741] 6.677465084194705e-05
[4.17408314 1.07094731] 6.0097184660515703e-05
[4.17408306 1.07094722] 5.408746531083866e-05
[4.17408298 1.07094714] 4.8678718040549386e-05
[4.17408292 1.07094707] 4.38108456636232e-05
[4.17408286 1.07094701] 3.942976062920491e-05
[4.1740828  1.07094695] 3.5486784190139235e-05
[4.17408275 1.07094689] 3.193810547693405e-05
[4.17408271 1.07094685] 2.8744294682388186e-05
[4.17408267 1.0709468 ] 2.5869865019915312e-05
[4.17408263 1.07094676] 2.328287837636998e-05
[4.1740826  1.07094673] 2.095459040508129e-05
[4.17408257 1.0709467 ] 1.885913123272889e-05
[4.17408255 1.07094667] 1.6973218042426106e-05
[4.17408253 1.07094664] 1.5275896164635222e-05
[4.1740825  1.07094662] 1.3748306463118147e-05
[4.17408249 1.0709466 ] 1.2373475779640278e-05
[4.17408247 1.07094658] 1.1136128161507551e-05
[4.17408245 1.07094657] 1.0022515325621178e-05
[4.17408244 1.07094655] 9.020263763958213e-06
[4.17408243 1.07094654] 8.118237382843987e-06
[4.17408242 1.07094653] 7.3064136099740275e-06
[4.17408241 1.07094651] 6.575772225897955e-06
[4.1740824  1.07094651] 5.9181950009792285e-06
[4.17408239 1.0709465 ] 5.326375504298625e-06
[4.17408238 1.07094649] 4.793737944882408e-06
[4.17408238 1.07094648] 4.314364149786832e-06
[4.17408237 1.07094647] 3.8829277264432105e-06
[4.17408236 1.07094647] 3.494634963524326e-06
[4.17408236 1.07094646] 3.145171475405804e-06
[4.17408236 1.07094646] 2.8306543192357545e-06
[4.17408235 1.07094645] 2.5475888794322023e-06
[4.17408235 1.07094645] 2.292829988712573e-06
[4.17408234 1.07094645] 2.0635469947081172e-06
[4.17408234 1.07094644] 1.8571922945662136e-06
[4.17408234 1.07094644] 1.6714730720789095e-06
[4.17408234 1.07094644] 1.504325774079711e-06
[4.17408234 1.07094644] 1.3538931969681648e-06
[4.17408233 1.07094644] 1.2185038724153025e-06
[4.17408233 1.07094643] 1.0966534826478562e-06
[4.17408233 1.07094643] 9.869881493484445e-07
[4.17408233 1.07094643] 8.882893304249515e-07
[4.17408233 1.07094643] 7.994603956585794e-07
[4.17408233 1.07094643] 7.195143691891446e-07
[4.17408233 1.07094643] 6.475629280187211e-07
[4.17408232 1.07094643] 5.828066282775071e-07
[4.17408232 1.07094642] 5.245259652603972e-07
[4.17408232 1.07094642] 4.7207334773654855e-07
[4.17408232 1.07094642] 4.248660212312865e-07
[4.17408232 1.07094642] 3.8237941206654656e-07
[4.17408232 1.07094642] 3.4414146963107337e-07
[4.17408232 1.07094642] 3.0972729336626635e-07
[4.17408232 1.07094642] 2.7875455871748084e-07
[4.17408232 1.07094642] 2.5087909920936314e-07
[4.17408232 1.07094642] 2.2579121452770435e-07
[4.17408232 1.07094642] 2.0321209840954908e-07
[4.17408232 1.07094642] 1.828908870837847e-07
[4.17408232 1.07094642] 1.6460179274758132e-07
[4.17408232 1.07094642] 1.4814160937702221e-07
[4.17408232 1.07094642] 1.3332743107856137e-07
[4.17408232 1.07094642] 1.1999469229298465e-07
[4.17408232 1.07094642] 1.0799524113585899e-07
[4.17408232 1.07094642] 9.719572615526817e-08
[4.17408232 1.07094642] 8.747615444075035e-08
[4.17408232 1.07094642] 7.872854826950011e-08
[4.17408232 1.07094642] 7.085568209466179e-08
[4.17408232 1.07094642] 6.377011695471628e-08
[4.17408232 1.07094642] 5.7393127680194086e-08
[4.17408232 1.07094642] 5.1653815159016166e-08
[4.17408232 1.07094642] 4.6488422849780995e-08
[4.17408232 1.07094642] 4.183956814758585e-08
[4.17408232 1.07094642] 3.7655614340071293e-08
[4.17408232 1.07094642] 3.3890069510083446e-08
[4.17408232 1.07094642] 3.050106740484884e-08
[4.17408232 1.07094642] 2.745094665079787e-08
[4.17408232 1.07094642] 2.4705822714418927e-08
[4.17408232 1.07094642] 2.223526306388371e-08
[4.17408232 1.07094642] 2.0011736212072055e-08
[4.17408232 1.07094642] 1.8010564648706633e-08
[4.17408232 1.07094642] 1.6209519752781463e-08
[4.17408232 1.07094642] 1.4588581010919568e-08
[4.17408232 1.07094642] 1.3129737811226125e-08
[4.17408232 1.07094642] 1.181675023229428e-08
[4.17408232 1.07094642] 1.0635105241351827e-08
[4.17408232 1.07094642] 9.571590523212458e-09
[4.17408232 1.07094642] 8.614412709736792e-09
[4.17408232 1.07094642] 7.752969691510336e-09
[4.17408232 1.07094642] 6.977691664728855e-09
[4.17408232 1.07094642] 6.279927551619311e-09
[4.17408232 1.07094642] 5.65195041535497e-09
[4.17408232 1.07094642] 5.086753966952846e-09
[4.17408232 1.07094642] 4.578079362205442e-09
[4.17408232 1.07094642] 4.120274281860069e-09
[4.17408232 1.07094642] 3.7082427241760343e-09
[4.17408232 1.07094642] 3.337406285193158e-09
[4.17408232 1.07094642] 3.0036817430826743e-09
[4.17408232 1.07094642] 2.7033079756946473e-09
[4.17408232 1.07094642] 2.432963157479154e-09
[4.17408232 1.07094642] 2.189671799911082e-09
[4.17408232 1.07094642] 1.9707027623208493e-09
[4.17408232 1.07094642] 1.7736393456054598e-09
[4.17408232 1.07094642] 1.5962786741347351e-09
[4.17408232 1.07094642] 1.4366587717557547e-09
[4.17408232 1.07094642] 1.292987537771636e-09
[4.17408232 1.07094642] 1.163689255640791e-09
[4.17408232 1.07094642] 1.0473131406772617e-09
[4.17408232 1.07094642] 9.42573031526931e-10

#kovariancni matice  
V = -np.linalg.inv(Hesse(par,x))

print('----------------------------------------------------')    
print('poloha x0 = {0:6.3f} +/- {1:6.3f} '.format(par[0],np.sqrt(V[0,0])))
print('vzdálenost l = {0:6.3f} +/- {1:6.3f}'.format(par[1],np.sqrt(V[1,1])))
print('korelace = {0:5.3f}'.format(V[0,1]/(np.sqrt(V[0,0]*V[1,1]))))
print('----------------------------------------------------')

----------------------------------------------------
poloha x0 =  4.174 +/-  0.148 
vzdálenost l =  1.071 +/-  0.155
korelace = -0.023
----------------------------------------------------

Newton-Raphsonův algoritmus¶