Rayleigh Statistics and GauCh for non-Gaussian Noise Analysis

This tutorial demonstrates the implementation of a comprehensive non-Gaussian noise analysis toolkit in GWexpy, based on research by Yamamura (2024) and Yamamoto (2015-2016).

Contents

Noise Simulation: Generating non-stationary and scattered light noise.
Rayleigh Statistics: Testing deviations from Rayleigh distribution.
GauCh (Modified KS Test): Advanced non-Gaussianity detection.
Student-t Indicator: Measuring distribution tail thickness.
Data Quality Flags: Generating veto segments automatically.
Evaluation & Visualization: ROC curves and composite dashboards.

[ ]:

import numpy as np
import matplotlib.pyplot as plt
from gwexpy.timeseries import TimeSeries
from gwexpy.noise import transient_gaussian_noise, scatter_light_noise, apply_line_mask
from gwexpy.statistics import compute_gauch, rayleigh_pvalue, compute_student_t_nu, to_segments
from gwexpy.plot.gauch_dashboard import plot_gauch_dashboard
%matplotlib inline

1. Noise Simulation

We simulate two types of non-Gaussian noise models used in KAGRA characterization.

Model I: Transient Gaussian noise (glitches).
Model II: Scattered light noise (stationary non-Gaussianity).

[ ]:

duration = 64.0
fs = 1024.0

# Model I: Glitch injection
ts_glitch = transient_gaussian_noise(duration, fs, A1=20.0, name='Glitchy Data')

# Model II: Scattered light
ts_scatter = scatter_light_noise(duration, fs, A2=1e-12, name='Scattered Light')

ts_glitch.plot()

2. Rayleigh Spectrogram

The Rayleigh statistic \(R\) measures the consistency of the ASD distribution with a Rayleigh distribution. We use the rayleigh_test method to get a p-value map.

[ ]:

rs_p = ts_glitch.rayleigh_test(fftlength=0.4, stride=0.5, overlap=0.2, n_samples=20)
rs_p.plot(yscale='log', title='Rayleigh Test p-value Map')
plt.colorbar(mappable=plt.gca().get_images()[-1] if plt.gca().get_images() else plt.gca().collections[-1], label='p-value')

3. GauCh (Modified KS Test)

GauCh is a more sensitive test for non-Gaussianity using a modified Kolmogorov-Smirnov test.

[ ]:

gauch_res = ts_glitch.gauch(fftlength=1.0, window=20)
gauch_res.pvalue_map.plot(yscale='log', title='GauCh p-value Map')
plt.colorbar(mappable=plt.gca().get_images()[-1] if plt.gca().get_images() else plt.gca().collections[-1], label='p-value')

4. Student-t Indicator

The Student-t indicator fits a Student-t distribution to FFT components and outputs the degree of freedom \(\nu\).

[ ]:

nu_spec = ts_glitch.student_t_spectrogram(fftlength=1.0, window=20)
nu_spec.plot(yscale='log', title='Student-t nu Map')
plt.colorbar(mappable=plt.gca().get_images()[-1] if plt.gca().get_images() else plt.gca().collections[-1], label='nu')

5. Data Quality Flags

We can automatically generate veto segments where the p-value is below a threshold.

[ ]:

dq_flag = to_segments(gauch_res.pvalue_map, alpha=0.001)
print(dq_flag)
dq_flag.plot()

6. Composite Dashboard

Finally, we can visualize everything in a single dashboard.

[ ]:

fig = plot_gauch_dashboard(ts_glitch, gauch_res)
plt.show()