Written by Rose A. Finn, updated on January 5, 2015
PURPOSE:
This program calculates the biweight center and scale for the LCS clusters
using existing programs from the astropy.stats package. Errors on the center
and scale are calculating using bootstrap resampling (1000 samples is the default).
CALLING SEQUENCE
from within ipython
% run ~/svnRepository/pythonCode/LCSbiweight.py
getbiweightall()
to see the results plotted with velocity histogram for one cluster
mkw11.plotvhist()
to see a multipanel plot for all clusters, type
plotall()
INPUT PARAMETERS
none
OUTPUT PARAMETERS
none
EXAMPLES
see calling sequence
PROCEDURE
The NSA catalog is cut to include only those galaxies with velocities within
4000 km/s of the central velocity (from the literature) and within a projected
radius of 1 degree. The biweight location and scale are calculated, and the
location is used as the median, and the biweight location and scale are recalculated.
This is repeated until the scale changes by less than 1 km/s. This typically requires
2 iterations.
REQUIRED PYTHON MODULES
numpy
pylab
astropy
scipy
In [1]:
from matplotlib import pyplot as plt
import numpy as np
import scipy.stats
from astropy.stats import biweight_location, biweight_midvariance, sigma_clip, bootstrap
from astropy import cosmology
from astropy.io import fits
import glob
%matplotlib inline
In [2]:
def centralbi(x, M=None):
x=np.array(x,'f')
if M is None:
M=np.median(x)
MAD=np.median(abs(x-M))
ui=((x-M)/(6*MAD))
top=np.sum((x-M)*((1-ui**2)**2))
bottom=np.sum((1-ui**2)**2)
cbi=M + (top/bottom)
#print self.clustername
# print(cbi)
#finds the biweight scale
n=len(x)
usbi=((x-M)/(9*MAD))
upper= np.sum(((x-M)**2)*((1-usbi**2)**4))
lower=np.sum((1-usbi**2)*(1-5*usbi**2))
sbi=np.sqrt(n)*((sqrt(upper))/(abs(lower)))
return cbi,sbi
def getbiweight(z):
scale_cut=3.
biweightscale=biweight_midvariance(z)
biweightlocation=biweight_location(z)
flag=np.abs(z-biweightlocation)/biweightscale < scale_cut
#flag=np.ones(len(z),'bool')
repeatflag=1
nloop=0
#print biweightlocation, biweightscale
oldbiweightscale=biweightscale
while repeatflag:
newdata=z[flag]
biweightscale=biweight_midvariance(newdata, M=biweightlocation)
biweightlocation=biweight_location(newdata, M=biweightlocation)
oldflag = flag
#flag=abs(z-biweightlocation)/biweightscale < scale_cut
nloop += 1
#print nloop, biweightlocation, biweightscale, len(newdata), sum(flag)
#if sum(flag == oldflag) == len(flag):
# repeatflag=0
# print 'nloop = ', nloop
if np.abs(biweightscale - oldbiweightscale) < 1.:
repeatflag=0
#print 'nloop = ', nloop
if nloop > 5:
repeatflag = 0
oldbiweightscale=biweightscale
#print nloop, biweightlocation, biweightscale
#flag=abs(z-biweightlocation)/biweightscale < 4.
#biweightscale=biweight_midvariance(z[flag],M=biweightlocation)
return biweightlocation, biweightscale
In [3]:
# read in sample.dat - has cluster, RA, Dec, recession velocity
names=[]
RA=[]
DEC=[]
vr=[]
infile1=open('sample.dat','r')
for line in infile1:
#print line
t=line.split()
names.append(t[0])
RA.append(t[1])
DEC.append(t[2])
vr.append(t[3])
infile1.close()
vr=np.array(vr,'f')
RA=np.array(RA,'f')
DEC=np.array(DEC,'f')
names=np.array(names)
print names
In [6]:
# read in fits AGC table
infiles=glob.glob('*_NSA.fits')
outfile=open('biweight_center_scale.dat','w')
outfile.write('# name center scale RA DEC \n')
for f in infiles:
clustername=f.split('_')[0]
clustervel=vr[names == clustername]
clusterRA=RA[names == clustername]
clusterDEC=DEC[names == clustername]
dat = fits.getdata(f)
# cut velocities to within +/- 4000 km/s
keepflag = (np.abs(dat.ZDIST*3.e5 - clustervel) < 4000.) & (np.sqrt((dat.RA-clusterRA)**2+(dat.DEC-clusterDEC)**2) < .75)
# calculate biweight center and scale
a,b=getbiweight(dat.ZDIST[keepflag]*3.e5)
print clustername,": center vel = %5.1f, scale = %5.1f, RA = %12.8f, DEC = %12.8f"%(a,b,clusterRA,clusterDEC)
outfile.write(clustername+" %5.1f %5.1f %12.8f %12.8f\n"%(a,b,clusterRA,clusterDEC))
outfile.close()
#print infiles
#dat.columns
In [4]:
In [4]: