import pyfstat import matplotlib.pyplot as plt import pandas as pd import numpy as np import os from tqdm import tqdm from oct2py import octave import glob from scipy.stats import rv_continuous, chi2 filenames = glob.glob("CollectedOutput/*.txt") plt.style.use('paper') Tspan = 100 * 86400 class maxtwoFinNoise_gen(rv_continuous): def _pdf(self, twoF, Ntrials): F = twoF/2.0 alpha = (1 + F)*np.exp(-F) a = Ntrials/2.0*F*np.exp(-F) b = (1 - alpha)**(Ntrials-1) return a*b df_list = [] for fn in filenames: df = pd.read_csv( fn, sep=' ', names=['dF0', 'dF1', 'twoF', 'runTime']) df['CLUSTER_ID'] = fn.split('_')[1] df_list.append(df) df = pd.concat(df_list) print 'Number of samples = ', len(df) fig, ax = plt.subplots() ax.hist(df.twoF, bins=50, histtype='step', color='k', normed=True, linewidth=1, label='Monte-Carlo histogram') maxtwoFinNoise = maxtwoFinNoise_gen(a=0) Ntrials_effective, loc, scale = maxtwoFinNoise.fit(df.twoF.values, floc=0, fscale=1) print 'Ntrials effective = {:1.2e}'.format(Ntrials_effective) twoFsmooth = np.linspace(0, df.twoF.max(), 1000) best_fit_pdf = maxtwoFinNoise.pdf(twoFsmooth, Ntrials_effective) ax.plot(twoFsmooth, best_fit_pdf, '-r', label=r'$p(2\mathcal{{F}}_{{\rm max}})$ for {} $N_{{\rm trials}}$' .format(pyfstat.texify_float(Ntrials_effective, d=2))) pval = 1e-6 twoFsmooth_HD = np.linspace( twoFsmooth[np.argmax(best_fit_pdf)], df.twoF.max(), 100000) best_fit_pdf_HD = maxtwoFinNoise.pdf(twoFsmooth_HD, Ntrials_effective) spacing = twoFsmooth_HD[1]-twoFsmooth_HD[0] print twoFsmooth_HD[np.argmin(np.abs(best_fit_pdf_HD - pval))], spacing ax.set_xlabel('$\widetilde{2\mathcal{F}}$') ax.set_xlim(0, 60) ax.legend(frameon=False, fontsize=6, loc=2) fig.tight_layout() fig.savefig('allsky_noise_twoF_histogram.png') from latex_macro_generator import write_to_macro write_to_macro('AllSkyMCNoiseOnlyMaximum', '{:1.1f}'.format(np.max(df.twoF)), '../macros.tex') write_to_macro('AllSkyMCNoiseN', len(df), '../macros.tex')