Source code for lsst.sims.maf.metrics.periodicDetectMetric

import numpy as np
from .baseMetric import BaseMetric
from lsst.sims.maf.utils import m52snr
import lsst.sims.utils as utils
import scipy

__all__ = ['PeriodicDetectMetric']


[docs]class PeriodicDetectMetric(BaseMetric): """Determine if we would be able to classify an object as periodic/non-uniform, using an F-test The idea here is that if a periodic source is aliased, it will be indistinguishable from a constant source, so we can find a best-fit constant, and if the reduced chi-squared is ~1, we know we are aliased. Parameters ---------- period : float (2) or array The period of the star (days). Can be a single value, or an array. If an array, amplitude and starMag should be arrays of equal length. amplitude : floar (0.1) The amplitude of the stellar variablility (mags). starMag : float (20.) The mean magnitude of the star in r (mags). sig_level : float (0.05) The value to use to compare to the p-value when deciding if we can reject the null hypothesis. SedTemplate : str ('F') The stellar SED template to use to generate realistic colors (default is an F star, so RR Lyrae-like) Returns ------- 1 if we would detect star is variable, 0 if it is well-fit by a constant value. If using arrays to test multiple period-amplitude-mag combinations, will be the sum of the number of detected stars. """ def __init__(self, mjdCol='observationStartMJD', periods=2., amplitudes=0.1, m5Col='fiveSigmaDepth', metricName='PeriodicDetectMetric', filterCol='filter', starMags=20, sig_level=0.05, SedTemplate='F', **kwargs): self.mjdCol = mjdCol self.m5Col = m5Col self.filterCol = filterCol if np.size(periods) == 1: self.periods = [periods] # Using the same magnitude for all filters. Could expand to fit the mean in each filter. self.starMags = [starMags] self.amplitudes = [amplitudes] else: self.periods = periods self.starMags = starMags self.amplitudes = amplitudes self.sig_level = sig_level self.SedTemplate = SedTemplate super(PeriodicDetectMetric, self).__init__([mjdCol, m5Col, filterCol], metricName=metricName, units='N Detected (0, %i)' % np.size(periods), **kwargs)
[docs] def run(self, dataSlice, slicePoint=None): result = 0 n_pts = np.size(dataSlice[self.mjdCol]) n_filt = np.size(np.unique(dataSlice[self.filterCol])) # If we had a correct model with phase, amplitude, period, mean_mags, then chi_squared/DoF would be ~1 with 3+n_filt free parameters. # The mean is one free parameter p1 = n_filt p2 = 3.+n_filt chi_sq_2 = 1.*(n_pts-p2) u_filters = np.unique(dataSlice[self.filterCol]) if n_pts > p2: for period, starMag, amplitude in zip(self.periods, self.starMags, self.amplitudes): chi_sq_1 = 0 mags = utils.stellarMags(self.SedTemplate, rmag=starMag) for filtername in u_filters: in_filt = np.where(dataSlice[self.filterCol] == filtername)[0] lc = amplitude*np.sin(dataSlice[self.mjdCol][in_filt]*(np.pi*2)/period) + mags[filtername] snr = m52snr(lc, dataSlice[self.m5Col][in_filt]) delta_m = 2.5*np.log10(1.+1./snr) weights = 1./(delta_m**2) weighted_mean = np.sum(weights*lc)/np.sum(weights) chi_sq_1 += np.sum(((lc - weighted_mean)**2/delta_m**2)) # Yes, I'm fitting magnitudes rather than flux. At least I feel kinda bad about it. # F-test for nested models Regression problems: https://en.wikipedia.org/wiki/F-test f_numerator = (chi_sq_1 - chi_sq_2)/(p2-p1) f_denom = 1. # This is just reduced chi-squared for the more complicated model, so should be 1. f_val = f_numerator/f_denom # Has DoF (p2-p1, n-p2) # https://stackoverflow.com/questions/21494141/how-do-i-do-a-f-test-in-python/21503346 p_value = scipy.stats.f.sf(f_val, p2-p1, n_pts-p2) if np.isfinite(p_value): if p_value < self.sig_level: result += 1 return result