Florent Leclercq,
Imperial Centre for Inference and Cosmology, Imperial College London,
florent.leclercq@polytechnique.org
import numpy as np
import scipy.linalg
from matplotlib import pyplot as plt
np.random.seed(123456)
%matplotlib inline
N=1000
epsilon=0.000000000001
t=np.array([100./(f*f) for f in range(1,N//2)])
signalcovar=np.concatenate([np.zeros(1),t,t[::-1],np.zeros(1)])
plt.plot(np.arange(N),signalcovar)
plt.show()
sqrtsignalcovar=np.diagflat(np.sqrt(signalcovar))
sqrtsignalcovarPix=np.fft.ifft(np.fft.fft(sqrtsignalcovar).T).T
plt.imshow(sqrtsignalcovarPix.real, cmap='RdBu')
plt.title("Signal covariance matrix")
plt.show()
mask=np.ones(N)
mask[400:800]=0.
plt.ylim(0,1.1)
plt.plot(np.arange(N),mask)
plt.fill_between([400,800],0.,1.1,facecolor='grey',alpha=0.3, linewidth=0.)
plt.title("Mask")
plt.show()
noisepower=0.00001
noisecovar=noisepower*np.concatenate([np.ones(N//2),10000*np.ones(N//2)])
plt.plot(np.arange(N),noisecovar)
plt.show()
invnoisecovarmat=np.diagflat(mask/noisecovar)
plt.imshow(np.log(invnoisecovarmat+epsilon), cmap='viridis')
plt.title("Inverse of the noise covariance matrix")
plt.show()
Data model: $d=s+n$
#The truth
normalsim=np.random.normal(0.,1.,N)
s=sqrtsignalcovarPix.real.dot(normalsim)
plt.xlim(0,N)
plt.ylim(-2,3)
plt.scatter(np.arange(N),s,marker='.',color='black',label="true signal")
plt.fill_between([400,800],-2,3,facecolor='grey',alpha=0.3, linewidth=0.)
plt.title("Signal")
plt.legend(loc='best')
plt.show()
#The noise, with infinite variance in masked regions
normalsim=np.random.normal(0.,1.,N)
n=normalsim/np.sqrt(np.diag(invnoisecovarmat)+epsilon)
d=s+n
plt.xlim(0,N)
plt.ylim(-2,3)
plt.scatter(np.arange(N),mask*n,marker='.',color='green',label="noise")
plt.scatter(np.arange(N),mask*d,marker='.',color='blue',label="masked data")
plt.fill_between([400,800],-2,3,facecolor='grey',alpha=0.3, linewidth=0.)
plt.title("Noise and data")
plt.legend(loc='best')
plt.show()
plt.xlim(0,N)
plt.ylim(-2,3)
plt.scatter(np.arange(N),s,marker='.',color='black',label="true signal")
plt.scatter(np.arange(N),mask*d,marker='.',color='blue',label="masked data")
plt.fill_between([400,800],-2,3,facecolor='grey',alpha=0.3, linewidth=0.)
plt.title("Signal and data")
plt.legend(loc="best")
plt.show()
Covariance of the Wiener Filter: \begin{equation} \mathrm{Cov}_\mathrm{WF} = (N^{-1}+S^{-1})^{-1} = S^{1/2}(I+S^{1/2}N^{-1}S^{1/2})^{-1}S^{1/2} \end{equation}
M=np.identity(N)+sqrtsignalcovarPix.dot(invnoisecovarmat).dot(sqrtsignalcovarPix)
CovWF=sqrtsignalcovarPix.dot(np.linalg.inv(M)).dot(sqrtsignalcovarPix)
CovWF=(CovWF+CovWF.T)/2
sqrtCovWF=scipy.linalg.sqrtm(CovWF)
plt.imshow(CovWF.real, cmap='RdBu')
plt.title("Covariance matrix of the Wiener Filter")
plt.show()
Mean of the Wiener posterior: \begin{equation} s_\mathrm{WF} = \mathrm{Cov}_\mathrm{WF} N^{-1} d \end{equation}
sWF=CovWF.dot(invnoisecovarmat).dot(d).real
plt.xlim(0,N)
plt.ylim(-2,3)
plt.scatter(np.arange(N),s,marker='.',color='black',label="true signal")
plt.scatter(np.arange(N),sWF,marker='.',color='red',label="reconstructed signal")
plt.fill_between([400,800],-2,3,facecolor='grey',alpha=0.3, linewidth=0.)
plt.title("Signal reconstruction")
plt.legend(loc='best')
plt.show()
Samples of the Wiener posterior: \begin{equation} s=s_\mathrm{WF}+\sqrt{C_\mathrm{WF}} \, G(0,1) \end{equation} so that $\left\langle s \right\rangle = s_\mathrm{WF}$ and $\mathrm{Cov}(s) = C_\mathrm{WF}$
cr1=sqrtCovWF.dot(np.random.normal(0.,1.,N)).real+sWF
cr2=sqrtCovWF.dot(np.random.normal(0.,1.,N)).real+sWF
cr3=sqrtCovWF.dot(np.random.normal(0.,1.,N)).real+sWF
plt.xlim(0,N)
plt.ylim(-2,3)
plt.scatter(np.arange(N),s,marker='.',color='black',label="true signal")
plt.scatter(np.arange(N),cr1,marker='.',color='yellow')
plt.scatter(np.arange(N),cr2,marker='.',color='pink')
plt.scatter(np.arange(N),cr3,marker='.',color='purple',label="constrained realization")
plt.fill_between([400,800],-2,3,facecolor='grey',alpha=0.3, linewidth=0.)
plt.title("Constrained realizations")
plt.legend(loc='best')
plt.show()
plt.xlim(0,N)
plt.ylim(-2,3)
plt.scatter(np.arange(N),s-cr1,marker='.',color='yellow')
plt.scatter(np.arange(N),s-cr2,marker='.',color='pink')
plt.scatter(np.arange(N),s-cr3,marker='.',color='purple')
plt.plot([0,N],[0,0],color='black',linestyle='--',linewidth=2)
plt.fill_between([400,800],-2,3,facecolor='grey',alpha=0.3, linewidth=0.)
plt.title("Residuals")
plt.show()