Gaussian random fields and local non-Gaussianity

Florent Leclercq,
Institute of Cosmology and Gravitation, University of Portsmouth,
[email protected]

In [1]:
import numpy as np
import scipy.linalg
from scipy.stats import norm
from matplotlib import pyplot as plt
%matplotlib inline
In [2]:
N=1000
epsilon=0.00000001

Generate white noise

In [3]:
rv=norm(loc=0.,scale=1.)
x=np.linspace(-5,5,100)
whitenoise=np.random.normal(0.,1.,N)
In [4]:
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', normed=True, color='black', histtype='step', lw=1.5)
plt.show()

Setup covariance and generate Gaussian signal

Covariance demo 1

In [5]:
covar1=np.array([[np.exp(-np.abs((i-j)/20.)) for i in xrange(N)] for j in xrange(N)])
In [6]:
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,5))
ax1.matshow(covar1, cmap='RdBu')
ax2.set_xlim(0,N)
ax2.plot(np.arange(N),covar1[N/2])
plt.show()
In [7]:
rootcovar1=scipy.linalg.sqrtm(covar1)
signal1=rootcovar1.dot(whitenoise).real
In [8]:
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', normed=True, histtype='step', lw=1.5)
plt.show()

Covariance demo 2

In [9]:
covar2=np.array([[1./(1+np.abs((i-j)/2.))**2 for i in xrange(N)] for j in xrange(N)])
In [10]:
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,5))
ax1.matshow(np.log(covar2+epsilon), cmap='RdBu')
ax2.set_xlim(0,N)
ax2.plot(np.arange(N),covar2[N/2])
plt.show()
In [11]:
rootcovar2=scipy.linalg.sqrtm(covar2)
signal2=rootcovar2.dot(whitenoise).real
In [12]:
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', normed=True, histtype='step', lw=1.5)
plt.show()

Covariance demo 3

In [13]:
covar3=np.array([[(1. if i<N/2 else 100.) if i==j else 0. for i in xrange(N)] for j in xrange(N)])
In [14]:
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,5))
ax1.matshow(np.log(covar3+epsilon), cmap='viridis')
ax2.set_xlim(0,N)
ax2.plot(np.arange(N),covar3[N/2])
plt.show()
In [15]:
rootcovar3=scipy.linalg.sqrtm(covar3)
signal3=rootcovar3.dot(whitenoise).real
In [16]:
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', normed=True, histtype='step', lw=1.5)
plt.show()

Covariance demo 4

In [17]:
covar4=np.array([[np.cos((i-j)/8.) for i in xrange(N)] for j in xrange(N)])
In [18]:
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()
In [19]:
rootcovar4=scipy.linalg.sqrtm(covar4)
signal4=rootcovar4.dot(whitenoise).real
In [20]:
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', normed=True, histtype='step', lw=1.5)
plt.show()

Generate non-Gaussian signal

In [21]:
fNL=0.2
NGsignal=whitenoise+fNL*whitenoise**2
In [22]:
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', normed=True, histtype='step', color='purple', lw=1.5)
plt.show()