To find the marginal distribution for the flux of the star, we want to calculate $P\left(F\mid\{D_{x,y}\},I\right)$, where $I$ is any background information like the model we are using, etc. The model we use to fit the data will involve several parameters, so we put them all together in $$ \vec{z}_0 = \left(\sigma_0,x_0,y_0,F,B\right) $$
We will describe these parameters later when they are calculated. The marginal distribution of the flux, $F$, is then
$$ P\left(F\mid\{D_{x,y}\},I\right) = \int p\left(\vec{z}_0 \mid \{D_{x,y}\}, I\right)d\sigma_0dx_0dy_0dB $$From Baye's theorem, assuming a uniform prior $p\left(\vec{z}_0\mid I\right)$, we have $$ p\left(\vec{z}_0 \mid \{D_{x,y}\}, I\right) = \frac{p\left(\{D_{x,y}\}\mid \vec{z}_0, I\right)}{\int p\left(\{D_{x,y}\}\mid \vec{z}_0, I\right) d^5\vec{z}_0} $$
Before we carry out the integral, we need to find the best fit parameters $\vec{z}_0$ given the data and the covariance matrix, which will give an estimate of the errors in the fit. To do this, we assume that at each pixel, the photon-counting distribution is Poisson. Moreover, since the counts are large, the Poisson distribution can be approximated as a Gaussian distribution. Specifically, at each pixel we have $$ p\left(\vec{z}_0 \mid \{D_{x,y}\}, I\right) \propto \exp\left[-\frac{(D_{x,y} - f(x,y,\vec{z}_0))^2}{2\sigma_{x,y}^2}\right] $$
We also assume the stellar point-spread function is azimuthally symmetric Gaussian, which means we fit the data to a function $f$ of the form $$ f(x,y,\vec{z}_0) = F\exp\left[-\frac{(x-x_0)^2 + (y-y_0)^2}{2\sigma_0}\right] + B $$
Using the least-squares fitting routine in scipy.optimize called curvefit, we find the best fit parameters and covariance matrix $\Sigma$ below. To do the fit, we assume $\sigma{x,y}^2 \propto D_{x,y}$. The results are shown below.
In [ ]:
In [2]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize as opt
from astropy.io import fits
In [3]:
#Import file and plot data
starfile = fits.open('hw3prob1-data.fits')
stardata = starfile[0].data
plt.contourf(stardata)
plt.colorbar()
Out[3]:
In [4]:
def gaussian_f(z,sig,x0,y0,F,B):
#Model function to fit the data
return F*np.exp(-1/(2*sig**2)*((z[0]-x0)**2+(z[1]-y0)**2)) + B
In [5]:
#Using curve_fit to fit stardata to gaussian_f. Initial values are chosen by inspection of the data.
#The parameters are in the list popt and the covariance matrix is in the matrix pcov.
x_vals = np.linspace(0,stardata.shape[0]-1,stardata.shape[0])
y_vals = np.linspace(0,stardata.shape[1]-1,stardata.shape[1])
x_vals,y_vals = np.meshgrid(x_vals,y_vals)
initial = (25.,125.,125.,80.,100.)
sigma_0 = np.sqrt(stardata)
popt, pcov = opt.curve_fit(gaussian_f, (x_vals.ravel(),y_vals.ravel()), stardata.ravel(), p0=initial
, sigma=sigma_0.ravel())
In [6]:
print(popt)
print('')
print(pcov)
In [7]:
#Plotting fit over data and zoomed in. Looks good!
plt.contourf(stardata)
plt.colorbar()
plt.contour(x_vals,y_vals,gaussian_f((x_vals,y_vals),popt[0],popt[1],popt[2],popt[3],popt[4]), colors='black')
plt.axis([100,150,100,160])
plt.xlabel('x [pix]')
plt.ylabel('y [pix]')
Out[7]:
Now that we have the best fit parameters and the covariance matrix, we can use fact that integrating the 5-dimensional Gaussian reduces to another Gaussian. The best fit parameters give us the means and the entries of the covariance matrix give us the variances. In particular, the marginal distribution for the flux of the star is $$ P\left(F\mid \{D_{x,y}\},I\right) = \frac{1}{\sqrt{2\pi\sigma_F^2}}\exp\left[-\frac{(F-F_0)^2}{2\sigma_F^2}\right] $$ where $F_0 = 26.19$ and $\sigma_F^2 = 0.64$, which are read off from the 4th entry of popt and the (4,4) entry of pcov. This is DN units. To convert to photoelectrons, we multiply both $F$ and $\sigma_F$ by 3.0 $e^-$/DN.
To get the marginal distribution for the position of the star, we use the same argument to say that the result of the marginalization over the parameters $(\sigma_0,F,B)$ is a 2-dimensional Gaussian. The result is $$ P\left(x,y\mid \{D_{x,y}\},I\right) = \frac{1}{2\pi\det\Sigma^\prime}\exp\left[-\frac{1}{2}X^T\Sigma^{\prime -1}X\right] $$ where $$ X = \left( \begin{array}{ccc} x - x_0\\ y - y_0 \end{array} \right) $$ , $x_0=125.97$, and $y_0= 132.45$. The matrix $\Sigma^\prime$ is given by the $$ \Sigma^\prime = \left( \begin{array}{ccc} 4.48\times 10^{-2} & 4.13\times 10^{-4}\\ 4.13\times 10^{-4} & 4.48\times 10^{-2} \end{array} \right) $$ These values are in units of pixels.
In [8]:
cov = np.array([pcov[1,1:3],pcov[2,1:3]])
print(cov)
In [9]:
P_array = np.array([[1/(2*np.pi*np.sqrt(np.linalg.det(cov)))*np.exp(-0.5*float(np.dot(np.dot(np.transpose(np.array([[xx-popt[1]],
[yy-popt[2]]])),
np.linalg.inv(cov)),np.array([[xx-popt[1]],
[yy-popt[2]]]))))
for xx in range(0,len(stardata))] for yy in range(0,len(stardata))])
plt.contourf(stardata)
plt.colorbar()
plt.contour(range(0,len(stardata)),range(0,len(stardata)), P_array, colors='black')
plt.axis([122,130,130,135])
plt.xlabel('x [pix]')
plt.ylabel('y [pix]')
Out[9]:
In [12]:
print(popt[3]*(popt[0]**2)*2*np.pi)
In [ ]: