estimate the std error from hessian in psffit

maoppy/example/fit_muse_simulated_psf.py

@@ -22,7 +22,7 @@ from maoppy.psffit import psffit
 from maoppy.instrument import muse_nfm
 
 npix = 128 # pixel size of PSF
-wvl = 600*1e-9 # wavelength [m]
+wvl = 800*1e-9 # wavelength [m]
 
 flux = 1e4
 ron = muse_nfm.ron
@@ -33,25 +33,28 @@ samp = muse_nfm.samp(wvl) # sampling (2.0 for Shannon-Nyquist)
 Pmodel = Psfao((npix,npix), system=muse_nfm, samp=samp)
 
 #%% Generate a MUSE-NFM PSF
-r0 = 0.15 # Fried parameter [m]
+r0 = 0.12 # Fried parameter [m]
 b = 1e-2 # Phase PSD background [rad² m²]
 amp = 3.0 # Phase PSD Moffat amplitude [rad²]
 alpha = 0.1 # Phase PSD Moffat alpha [1/m]
-ratio = 1.0
-theta = 0.0
+ratio = 1.2 # major/minor axis
+theta = np.pi/3 # angle [rad]
 beta = 1.6 # Phase PSD Moffat beta power law
 
 param = [r0,b,amp,alpha,ratio,theta,beta]
 
-psf = Pmodel(param,dx=0,dy=0)
+dxtrue = 0
+dytrue = 0
+
+psf = Pmodel(param,dx=dxtrue,dy=dytrue)
 
 noise = ron*np.random.randn(npix,npix)
 image = flux*psf + bckgd + noise
 
 #%% Fit the PSF with Psfao
-guess = [0.145,2e-2,1.2,0.08,ratio,theta,1.5]
+guess = [0.145,2e-2,1.2,alpha,1,0,1.5]
 w = np.ones_like(image)/ron**2.0
-fixed = [False,False,False,False,True,True,False]
+fixed = [False,False,False,True,False,False,False]
 
 out = psffit(image, Pmodel, guess, weights=w, fixed=fixed)
 
@@ -59,14 +62,20 @@ flux_fit, bck_fit = out.flux_bck
 fitao = flux_fit*out.psf + bck_fit
 
 #%% Plots
-print(31*'-')
-print("%9s %10s %10s"%('','MUSE-NFM','Psfao'))
-print(31*'-')
-print("%9s %10.4g %10.4g"%('Flux [e-]',flux,flux_fit))
-print("%9s %10.4g %10.4g"%('Bck [e-]',bckgd,bck_fit))
-print("%9s %10.2f %10.2f"%('r0 [cm]',r0*100,out.x[0]*100))
-print("%9s %10.2f %10.2f"%('A [rad²]',amp,out.x[2]))
-print(31*'-')
+print(42*'-')
+print("%9s %10s %10s %10s"%('','TRUE','ESTIM','STD'))
+print(42*'-')
+for i in range(len(param)):
+    if not fixed[i]:
+        print("%9s %10.3f %10.3f %10.3f"%(Pmodel.param_name[i],param[i],out.x[i],out.x_std[i]))
+# print("%9s %10.2f %10.2f %10.2f"%('r0 [cm]',r0*100,out.x[0]*100,out.x_std[0]*100))
+# print("%9s %10.2f %10.2f %10.2f"%('A [rad²]',amp,out.x[2],out.x_std[2]))
+print(' -'*21)
+print("%9s %10.3g %10.3g"%('Flux [e-]',flux,flux_fit))
+print("%9s %10.3g %10.3g"%('Bck [e-]',bckgd,bck_fit))
+print("%9s %10.2f %10.2f %10.2f"%('dx [pix]',dxtrue,out.dxdy[0],out.dxdy_std[0]))
+print("%9s %10.2f %10.2f %10.2f"%('dy [pix]',dytrue,out.dxdy[1],out.dxdy_std[1]))
+print(42*'-')
 
 vmin,vmax = np.arcsinh([image.min(),image.max()])
 

maoppy/psffit.py

@@ -135,7 +135,7 @@ def psffit(psf,model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
             x, dxdy = mini2input(y)
             mm = model(x, dx=dxdy[0], dy=dxdy[1], **kwargs)
             amp, bck = lsq_flux_bck(mm, psf, weights, background=flux_bck[1], positive_bck=positive_bck)
-            return np.reshape(sqw * (amp*mm + bck - psf), np.size(psf))
+            return 0.5*np.reshape(sqw * (amp*mm + bck - psf), np.size(psf))
     
     cost = CostClass()
     
@@ -178,6 +178,11 @@ def psffit(psf,model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
     print("") #finish line of "-"
     
     result.x, result.dxdy = mini2input(result.x) #split output between x and dxdy
+    
+    std = 1/np.sqrt(np.diag(result.jac.T @ result.jac))
+    result.x_std, result.dxdy_std = mini2input(std)
+    if fixed is not None:
+        result.x_std[fixed] = 0
 
     m = model(result.x,dx=result.dxdy[0],dy=result.dxdy[1])
     amp, bck = lsq_flux_bck(m, psf, weights, background=flux_bck[1], positive_bck=positive_bck)