<jupyter_start><jupyter_text>Study of Variable StarsIn this notebook, we study the variability of the star _RR Lyrae_. We first download the data (flux, error in flux and magnitude) of the object from the Gaia data archive.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.optimize import minimize as mz
from scipy.optimize import curve_fit as cft
from scipy.signal import lombscargle as lsg
import emcee
import corner<jupyter_output><empty_output><jupyter_text>Below, we read data from the `.csv` file. We use `pandas` to do so.<jupyter_code>data = pd.read_csv('anonymous1619179984702__1.csv', delimiter=',')
data
time, flx, flx_err, mag = data['#time'].to_numpy(), data['#flux'].to_numpy(), data['#flux_over_error'].to_numpy(), data['#mag'].to_numpy()<jupyter_output><empty_output><jupyter_text>Note, that the data do not contain the error in magnitude, which we can compute from the error in the flux. We know that,$$M =-2.5\log{F} + C_0$$where, $M$, $F$ and $C_0$ corresponds to the magnitude, flux and the reference point, respectively. To find the error in $M$, we can derive above equation to find that,$$dM = 2.5\frac{1}{\ln{10}}\frac{dF}{F}$$$$\Rightarrow dM = 1.086\frac{dF}{F}$$below, we use this formula to find the error in magnitudes.  Variability visulizationBelow we present the plots of time v/s flux and time v/s magnitude to show the variability of the star.<jupyter_code>mag_err = 1.086*flx_err/flx
# Un-weighted mean, median: Flux
mn_fl, md_fl = np.average(flx), np.median(flx)
# Weighted mean: flux
mn_fl_wt = np.average(flx, weights=flx_err**-2)
#IQR
fl_iqr25, fl_iqr75 = np.quantile(flx, np.array([0.25,0.75]))

print('Un-weighted mean of the flux: ', mn_fl)
print('Weighted mean of the flux: ', mn_fl_wt)
print('Median of the flux: ', md_fl)
print('IQR: ' + '[' + str(fl_iqr25) + ', ' + str(fl_iqr75) + ']')
plt.figure(figsize=(16,9))
plt.errorbar(time, flx, yerr=flx_err, fmt='o', mfc='white', elinewidth=3, mew=2)
plt.plot(time, mn_fl*np.ones(len(flx)), label='Mean Flux')
plt.plot(time, mn_fl_wt*np.ones(len(flx)), label='Weighted Mean Flux')
plt.plot(time, md_fl*np.ones(len(flx)), label='Median Flux')
plt.fill_between(time, fl_iqr25*np.ones(len(flx)), fl_iqr75*np.ones(len(flx)), alpha=0.2, label='IQR')
plt.xlabel('Time')
plt.ylabel('Flux')
plt.legend(loc='best')
plt.grid()
# Un-weighted mean, median: Magnitude
mn_mag, md_mag = np.average(mag), np.median(mag)
# Weighted mean: flux
mn_mag_wt = np.average(mag, weights=mag_err**-2)
#IQR
mag_iqr25, mag_iqr75 = np.quantile(mag, np.array([0.25,0.75]))

print('Un-weighted mean of the flux: ', mn_mag)
print('Weighted mean of the flux: ', mn_mag_wt)
print('Median of the flux: ', md_mag)
print('IQR: ' + '[' + str(mag_iqr25) + ', ' + str(mag_iqr75) + ']')
plt.figure(figsize=(16,9))
plt.errorbar(time, mag, yerr=mag_err, fmt='o', mfc='white', elinewidth=3, mew=2)
plt.plot(time, mn_mag*np.ones(len(mag)), label='Mean Magnitude')
plt.plot(time, mn_mag_wt*np.ones(len(mag)), label='Weighted Mean Magnitude')
plt.plot(time, md_mag*np.ones(len(mag)), label='Median Magnitude')
plt.fill_between(time, mag_iqr25*np.ones(len(mag)), mag_iqr75*np.ones(len(mag)), alpha=0.2, label='IQR')
plt.xlabel('Time')
plt.ylabel('Magnitude')
plt.legend(loc='best')
plt.ylim([np.max(mag)+0.5, np.min(mag)-0.5])
plt.grid()<jupyter_output><empty_output><jupyter_text>Un-weighted mean or weighted mean?It is evident from both of above plots that magnitudes are much more regular -- we would work with magnitudes. However, it is not clear if we want to use the weighted mean or the un-weighted mean in our usage. To test that which is better to use, we perform a small simulation. We simulte an observation with true flux 50, and we take 15 observations -- since this process follows the Poisson distribution, we use `numpy.random.poisson` module to sample Poisson distribution of our observations. We then find the weighted mean and un-weighted mean of this observation. We repeat this procedure for many times and then we see which type of mean is more near (or more times occur) to the true mean value of flux.<jupyter_code>num_obs = 15
flux_true = 50

# Poisson distribution of flux to test the function
poi_flx = np.random.poisson(flux_true, num_obs)
print('Poisson distribution: ', poi_flx)
poi_flx_err = poi_flx**0.5
print('Error in it: ', poi_flx_err)
print('Un-weighted mean of the sampled distribution: ', np.average(poi_flx))
print('Weighted mean of the sampled distribution: ', np.average(poi_flx, weights=poi_flx_err**-2))
means = np.zeros(10000)
wt_means = np.zeros(10000)

for i in range(len(means)):
    poi_flxs = np.random.poisson(flux_true, num_obs)
    poi_flxs_err = poi_flxs**0.5
    means[i] = np.average(poi_flxs)
    wt_means[i] = np.average(poi_flxs, weights=poi_flxs_err**-2)
plt.figure(figsize=(16,9))
plt.hist(means, bins=20, histtype='step', linewidth=4, label='Un-weighted Mean')
plt.hist(wt_means, bins=20, histtype='step', color='g', linewidth=4, label='Weighted mean')
plt.axvline(50, c='r', label='True Flux', linewidth=4)
plt.xlabel('Flux', fontsize=16)
plt.legend(fontsize=16, loc='best')
plt.xlabel('Flux')
plt.ylabel('Conuts')
plt.grid()
plt.show()<jupyter_output><empty_output><jupyter_text>It is now clear from above figure that un-weighted mean better describe the true flux than the weighted mean.  Lomb-scargle analysisWe now try to find the periodicity in our time series data by using the Lomb-Scargle algorithm. Below, we present this analysis.<jupyter_code>days = np.linspace(1,365*5,10000)
freq = 2*np.pi/days
power = lsg(time, flx, freq)
plt.figure(figsize=(16,9))
plt.plot(days, power, c='g', linewidth=3)
plt.xlabel('Days')
plt.ylabel('Power')
plt.grid()<jupyter_output><empty_output><jupyter_text>Although there is a local maximum near $\sim 300 \ days$, it certainly not the global maximum!!!  Modelling using scipyBelow, we try to fit a sinusoidal function to the variablity of magnitude to find the period of the variablity of the star.<jupyter_code>def func(time, k, con, per, phi):
    ab = k*np.sin(2*np.pi*time/per + phi) + con
    return ab

xinit = np.array([1, mn_mag, 300, np.pi/2])

soln, pcov = cft(func, time, mag, sigma=mag_err, p0=xinit)
tt1 = np.linspace(np.min(time), np.max(time), 1000)

plt.figure(figsize=(16,9))
plt.errorbar(time, mag, yerr=mag_err, fmt='o', mfc='white', elinewidth=3, mew=2, label='Data')
plt.plot(tt1, func(tt1, *soln), c='orangered', lw=3, label='Best fitted model')
plt.xlabel('Days')
plt.ylabel('Magnitude')
plt.legend(loc='best')
plt.ylim([np.max(mag)+0.5, np.min(mag)-0.5])
plt.grid()
print('Period of the signal (in days): {:.4f} +/- {:.4f}'.format(soln[2], pcov[2][2]))<jupyter_output>Period of the signal (in days): 328.6282 +/- 11.2796
{'fraction_non_alphanumeric': 0.15078658355595131, 'fraction_numerical': 0.027753042445829622, 'mean_word_length': 5.904713114754099, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '166816', 'license': 'MIT', 'n_tokens_mistral': 2716, 'n_tokens_neox': 2500, 'n_words': 691, 'path': 'Large_scale_surveys/rr_lyrae.ipynb', 'repo_name': 'Jayshil/Data_analysis', 'size': '6738'}