#!/usr/bin/env python
# coding: utf-8
# # Gaussian random fields and local non-Gaussianity
# Florent Leclercq,
# Imperial Centre for Inference and Cosmology, Imperial College London,
# florent.leclercq@polytechnique.org
# In[1]:
import numpy as np
import scipy.linalg
from scipy.stats import norm
from matplotlib import pyplot as plt
np.random.seed(123456)
get_ipython().run_line_magic('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', density=True, color='black', histtype='step', lw=1.5)
plt.show()
# ## Setup covariance and generate Gaussian signal
# ### Covariance demo 1
# $C_{ij} = \exp\left(-\dfrac{|i-j|}{20}\right)$
# In[5]:
covar1=np.array([[np.exp(-np.abs((i-j)/20.)) for i in range(N)] for j in range(N)])
# In[6]:
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()
# 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', density=True, histtype='step', lw=1.5)
plt.show()
# ### Covariance demo 2
# $C_{ij}= \dfrac{1}{(1+|i-j|/2)^2}$
# In[9]:
covar2=np.array([[1./(1+np.abs((i-j)/2.))**2 for i in range(N)] for j in range(N)])
# In[10]:
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()
# 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', density=True, histtype='step', lw=1.5)
plt.show()
# ### Covariance demo 3
# $C_{ii}=\left\{
# \begin{array}{ll}
# 1~\mathrm{if}~i
# $C_{ij}=0$ for $i\neq j$
# In[13]:
covar3=np.array([[(1. if i