Florent Leclercq,
Imperial Centre for Inference and Cosmology, Imperial College London,
florent.leclercq@polytechnique.org
import numpy as np
import scipy.linalg
from scipy.stats import norm
from matplotlib import pyplot as plt
np.random.seed(123456)
%matplotlib inline
N=1000
epsilon=0.00000001
rv=norm(loc=0.,scale=1.)
x=np.linspace(-5,5,100)
whitenoise=np.random.normal(0.,1.,N)
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,5))
ax1.set_xlim(0,N)
ax1.scatter(np.arange(N),whitenoise,marker='.',color='black')
ax2.set_xlim(-5,5)
ax2.plot(x, rv.pdf(x), color='black', lw=2)
ax2.hist(whitenoise, bins='auto', density=True, color='black', histtype='step', lw=1.5)
plt.show()
$C_{ij} = \exp\left(-\dfrac{|i-j|}{20}\right)$
covar1=np.array([[np.exp(-np.abs((i-j)/20.)) for i in range(N)] for j in range(N)])
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,5))
ax1.imshow(covar1, cmap='RdBu')
ax2.set_xlim(0,N)
ax2.plot(np.arange(N),covar1[N//2])
plt.show()
rootcovar1=scipy.linalg.sqrtm(covar1)
signal1=rootcovar1.dot(whitenoise).real
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,5))
ax1.set_xlim(0,N)
ax1.scatter(np.arange(N),signal1,marker='.',color='blue')
ax2.set_xlim(-5,5)
ax2.plot(np.linspace(-5,5,100), norm(loc=0.,scale=1.).pdf(np.linspace(-5,5,100)), color='black', lw=2)
ax2.hist(signal1, bins='auto', density=True, histtype='step', lw=1.5)
plt.show()
$C_{ij}= \dfrac{1}{(1+|i-j|/2)^2}$
covar2=np.array([[1./(1+np.abs((i-j)/2.))**2 for i in range(N)] for j in range(N)])
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,5))
ax1.imshow(np.log(covar2+epsilon), cmap='RdBu')
ax2.set_xlim(0,N)
ax2.plot(np.arange(N),covar2[N//2])
plt.show()
rootcovar2=scipy.linalg.sqrtm(covar2)
signal2=rootcovar2.dot(whitenoise).real
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,5))
ax1.set_xlim(0,N)
ax1.scatter(np.arange(N),signal2,marker='.',color='blue')
ax2.set_xlim(-5,5)
ax2.plot(np.linspace(-5,5,1000), norm(loc=0.,scale=1.).pdf(np.linspace(-5,5,1000)), color='black', lw=2)
ax2.hist(signal2, bins='auto', density=True, histtype='step', lw=1.5)
plt.show()
$C_{ii}=\left\{
\begin{array}{ll}
1~\mathrm{if}~i<N/2\\
100~\mathrm{otherwise}
\end{array}
\right.$
$C_{ij}=0$ for $i\neq j$
covar3=np.array([[(1. if i<N/2 else 100.) if i==j else 0. for i in range(N)] for j in range(N)])
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,5))
ax1.imshow(np.log(covar3+epsilon), cmap='viridis')
ax2.set_xlim(0,N)
ax2.plot(np.arange(N),covar3[N//2])
plt.show()
rootcovar3=scipy.linalg.sqrtm(covar3)
signal3=rootcovar3.dot(whitenoise).real
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,5))
ax1.set_xlim(0,N)
ax1.scatter(np.arange(N),signal3,marker='.',color='blue')
ax2.set_xlim(-30,30)
ax2.plot(np.linspace(-5,5,1000), norm(loc=0.,scale=1.).pdf(np.linspace(-5,5,1000)), color='black', lw=2)
ax2.plot(np.linspace(-30,30,1000), norm(loc=0.,scale=10.).pdf(np.linspace(-30,30,1000)), color='black', lw=2)
ax2.hist(signal3, bins='auto', density=True, histtype='step', lw=1.5)
plt.show()
$C_{ij}=\cos\left(\dfrac{i-j}{8}\right)$
covar4=np.array([[np.cos((i-j)/8.) for i in range(N)] for j in range(N)])
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,5))
ax1.matshow(covar4, cmap='RdBu')
ax2.set_xlim(0,N)
ax2.plot(np.arange(N),covar4[N//2])
plt.show()
rootcovar4=scipy.linalg.sqrtm(covar4)
signal4=rootcovar4.dot(whitenoise).real
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,5))
ax1.set_xlim(0,N)
ax1.scatter(np.arange(N),signal4,marker='.',color='blue')
ax2.set_xlim(-5,5)
ax2.plot(np.linspace(-5,5,1000), norm(loc=0.,scale=1.).pdf(np.linspace(-5,5,1000)), color='black', lw=2)
ax2.hist(signal4, bins='auto', density=True, histtype='step', lw=1.5)
plt.show()
Histograms of Gaussian random fields are not necessarily Gaussian!
A non-Gaussian signal, $s = \Phi + f_\mathrm{NL} \Phi^2$ where $\Phi$ is a Gaussian random field. In cosmology, this is called "local-type" non-Gaussianity. The 1-point pdf is skewed.
fNL=0.2
NGsignal=whitenoise+fNL*whitenoise**2
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,5))
ax1.set_xlim(0,N)
ax1.scatter(np.arange(N),NGsignal,marker='.',color='purple')
ax2.set_xlim(-5,5)
ax2.plot(x, rv.pdf(x), color='black', lw=2)
ax2.hist(NGsignal, bins='auto', density=True, histtype='step', color='purple', lw=1.5)
plt.show()