import numpy as np
import matplotlib.pyplot as plt

hustota pravděpodobnosti Breit-Wignerova rozdělení

def breit_wigner(x,x0,l):
    return 1/np.pi*l/(l**2+(x-x0)**2)

Simulace dat majáku

x0=5  #poloha majáku
l=10  #vzdálenost majáku od pobřeží 
Ndata=10000 #počet simulovaných hodnot
data=l*np.tan(np.pi*np.random.random(Ndata))+x0
plt.hist(data,bins=100,range=(-100,100),density=True,label='simulace')
xp=np.arange(-100,100,0.1)
plt.plot(xp,breit_wigner(xp,x0,l),c='red',label='Breit-Wignerovo rozdělení')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.show()

Nyní vypočítáme hodnoty logaritmu věrohodnostní funkce na 2D mřížce v prostoru parametrů x0, l a najdeme kde v tomto prostoru parametrů je logaritmus věrohodnostní funkce maximální

x0f=np.linspace(-100,100,100)
lf=np.linspace(1,50,100)
l_mesh,x_mesh=np.meshgrid(lf,x0f) #vytvoří 2D mřížku v prostoru parameterů x0,l
logL=np.zeros([np.size(x0f),np.size(lf)])
for i in range(Ndata):
     logL+=np.log(breit_wigner(data[i],x_mesh,l_mesh))
imax,jmax=np.where(logL==np.max(logL))

Nakreslíme vrstevnicový graf logaritmu věrohodnostní funkce

plt.contour(x_mesh,l_mesh,logL,levels=100)
plt.scatter(x0f[imax],lf[jmax],c='red',marker='+')
plt.xlabel('poloha x0')
plt.ylabel('vzdálenost l')
plt.show()
print(f"poloha majáku  x0 = {x0f[imax]}, vzdálenost l =  {lf[jmax]}")

poloha majáku  x0 = [5.05050505], vzdálenost l =  [9.90909091]