Problém majáku

import numpy as np
import matplotlib.pyplot as plt

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

parametry pro simulaci

N=1000 #počet simulací
l_true=1 #vzdálenost majáku od pobřeží
x0_true=4 #x-ová souřadnice majáku

simulace

x=np.empty(N) #pole, kde budou souřadnice detektorů, které registrovaly signál majákuy
#simulace náhodné proměnné z Breit-Wignerova rozdělení metodou inverzní funkce
x=l_true*np.tan((np.random.random(N)-0.5)*np.pi)+x0_true
#x,y souřadnice lorenziánu pro vykreslení teoretické pdf
x_plot=np.linspace(x0_true-10*l_true,x0_true+10*l_true,1000)
y_plot=lorenzian(x_plot,x0_true,l_true)

fig,ax=plt.subplots()
ax.hist(x,bins=2*int(np.sqrt(N)),range=(x0_true-10*l_true,x0_true+10*l_true),density=True,label='simulace')
ax.plot(x_plot,y_plot,c='red',label='Breit-Wignerovo rozdělení')
ax.set_xlabel('x')
ax.set_ylabel('pdf')
ax.legend()
plt.show()

odhad parametrů

#síť pro likelihood funkci
Np=100 #počet bodů sítě v jednom směru
l=np.linspace(0.01,3*l_true,Np)
x0=np.linspace(x0_true-10*l_true,x0_true+10*l_true,Np)
x0_mesh,l_mesh=np.meshgrid(x0,l)
#print(x0_mesh)
Likelihood=np.zeros([Np,Np])
#výpočet log věrohodnostní funkce
for i in range(N):
    Likelihood+=np.log(lorenzian(x[i],x0_mesh,l_mesh))
#nalezení polohy maxima log věrohodnostní funkce 
jmax,imax=np.unravel_index(np.argmax(Likelihood),Likelihood.shape)
x0_fit=x0[imax]
l_fit=l[jmax]
print(f" odhad polohy majáku:      {x0_fit:6.2f}")
print(f" odhad vzdálenosti majáku: {l_fit:6.2f}")

 odhad polohy majáku:        4.03
 odhad vzdálenosti majáku:   0.93

cmap=plt.cm.jet
fig,ax=plt.subplots()
ax.pcolor(x0,l,Likelihood,cmap=cmap)
cnt=ax.contour(x0,l,Likelihood,levels=100,colors='k')
cnt.set_linestyle('solid')
cnt.set_linewidth(0.5)
ax.scatter(x0_true,l_true,s=10,c='black',label='skutečná hodnota')
ax.scatter(x0_fit,l_fit,s=100,marker='+',c='blue',label='odhad')
ax.set_xlabel('x0')
ax.set_ylabel('l')
ax.legend(loc='upper right')
plt.title("Věrohodnostní funkce")
plt.show()

odhad neurčitostí parametrů

#aproximace posteriorní pdf kvadratickou formou v okolí maxima
A=0
B=0
C=0
for i in range(N):
    A=A+(l_fit**2-(x[i]-x0_fit)**2)/(l_fit**2+(x[i]-x0_fit)**2)**2
    B=B+(l_fit**2-(x[i]-x0_fit)**2)/(l_fit**2+(x[i]-x0_fit)**2)**2
    C=C+(x[i]-x0_fit)/(l_fit**2+(x[i]-x0_fit)**2)**2
A=-2*A
B=-N/l_fit**2+2*B
C=-4*l_fit*C    
#print(A,B,C)
sigma_x0=np.sqrt(B/(C**2-A*B))
sigma_l=np.sqrt(A/(C**2-A*B))
cov=-C/(C**2-A*B)
rho=cov/(sigma_x0*sigma_l)
print('odhad polohy majáku     x0 = {0:5.2f} +/- {1:5.2f}'.format(x0_fit,sigma_x0))
print('odhad vzdálenosti majáku l = {0:5.2f} +/- {1:5.2f}'.format(l_fit,sigma_l))
print('korelace = {0:6.4f}'.format(rho))

odhad polohy majáku     x0 =  4.03 +/-  0.04
odhad vzdálenosti majáku l =  0.93 +/-  0.04
korelace = -0.0118

marginální posteriorní pdf

#marginal posterior pdfs
f_marginal_l=np.zeros(Np)
f_marginal_x0=np.zeros(Np)
sum_marginal_l=0
sum_marginal_x0=0
for j in range(Np):
    for k in range(Np):
        f_marginal_l[j]+=Likelihood[j,k]
        f_marginal_x0[j]+=Likelihood[k,j]
#normalizace
f_marginal_x0=f_marginal_x0/abs(np.sum(f_marginal_x0))
f_marginal_l=f_marginal_l/abs(np.sum(f_marginal_l))

fig,ax=plt.subplots()
ax.plot(x0,f_marginal_x0)
ax.scatter(x0_fit,np.max(f_marginal_x0),c='red',s=20)
ax.set_xlabel('x0',fontsize=12)
ax.set_ylabel('ln f(x0|Dn)',fontsize=12)
plt.title('marginální posteriorní pdf - souřadnice x0')
plt.show()

fig,ax=plt.subplots()
ax.plot(l,f_marginal_l)
ax.scatter(l_fit,np.max(f_marginal_l),c='red',s=20)
ax.set_xlabel('l',fontsize=12)
ax.set_ylabel('ln f(l|Dn)',fontsize=12)
plt.title('marginální posteriorní pdf - vzdálenost l')
plt.show()