python howtos

On this page... (hide)

  1.   1.  open a fits file and plot spectrum (loglog)
  2.   2.  read a likehood file and plot the points with errors
  3.   3.  covariance and correlation

1.  open a fits file and plot spectrum (loglog)

from astropy.io import fits
from matplotlib import pyplot as plt

h= fits.open("output.fits")
pkdata=h[3].data

k=pkdata['k']
pk=pkdata['pklin']

plt.loglog(k,pk)

2.  read a likehood file and plot the points with errors

The file is in plain text, columns corresponds to "k_i, Pk_i, sig(Pk_i)" 

import numpy as np
from matplotlib import pyplot as plt

k,pk,sigpk=np.loadtxt("pk_0.5.dat",skiprows=1,unpack=True)

plt.errorbar(k,pk,sigpk,fmt='.')

3.  covariance and correlation

When we have a set x_i,i=1,..,N of random variables the generalization of the variance of a parameter is given by the covariance matrix. 

 cov(x_i,x_j)=E[(x_i-\bar x_i)(x_j-\bar x_j)

where \bar x_i is the mean of x_i.

  • what are the diagonal terms?
  • which kind of matrix is it?

The correlation matrix is the covariance matrix divided by some diagonal elements.

cor(x_i,x_j)= \frac{cov(x_i,x_j)}{\sigma_i \sigma_j}

where \sigma_i is the standard deviation of variable x_i.

  • what are its diagonal terms?

By construction all the elements should lie in the range [-1,1]. This a convenient way to see how much pairs of variable are correlated: positive values= correlated, negative=anti-correlated, 0=uncorrelated.

To code it we could use a double index loop, here is a more pythonic way: 

def cov2cor(c):
    d=diag(c)
    return c/sqrt(d[:,None]*d[None,:])

back