150326 commit

%% Cell type:markdown id:91330533 tags:

**What is this?**

*This jupyter notebook is part of a collection of notebooks on various topics discussed during the Time Domain Astrophysics course delivered by Stefano Covino at the [Università dell'Insubria](https://www.uninsubria.eu/) in Como (Italy). Please direct questions and suggestions to [stefano.covino@inaf.it](mailto:stefano.covino@inaf.it).*

%% Cell type:markdown id:915ee876 tags:

**This is a `textual` notebook**

%% Cell type:markdown id:53194e25 tags:

![Time Domain Astrophysics](Pics/TimeDomainBanner.jpg)

%% Cell type:markdown id:3a9502d7 tags:

# Baeysian view of the LS Periodogram

***

- What we want to be able to do is to detect variability and measure the period in the face of both noisy and incomplete data. Instead we'll use Fourier decomposition to get a more useful tool for actual data analysis.

- For a periodic signal we have:

$$y(t+P)=y(t),$$ where $P$ is the period.

- We can create a *phased light curve* that plots the data as function of phase:

$$\phi=\frac{t}{P} − {\rm int}\left(\frac{t}{P}\right),$$

- where ${\rm int}(x)$ returns the integer part of $x$.

%% Cell type:markdown id:6deddd35-e973-43a0-a658-8a2475af688a tags:

### A Single Sinusoid

***

- Let's take the case where the data are drawn from a single sinusoidal signal:

$$y(t)=A \sin(\omega t+\phi)+\epsilon$$

- and determine whether or not the data are indeed consistent with periodic variability and, if so, what is the period.

- This model is **non-linear** in the frequency term, $\omega$ and the phase, $\phi$ and therefore We rewrite the argument as $\omega(t−t_0)$ (reexpressing the phase term) and use trigonometrics identies to rewrite the model as:

$$y(t)=a \sin(\omega t)+b \cos(\omega t)$$

- where

$$A=(a^2+b^2)^{1/2} \text{ and } \phi=\tan^{−1}(b/a)$$

- The model is now linear with respect to coefficients $a$ and $b$ (and nonlinear only with respect to frequency, $\omega$).

%% Cell type:markdown id:270f9483-19e5-44d0-a873-c32998db252e tags:

- Assuming constant uncertainties on the data, we can write a likelihood function down:

$$L =\prod^N_{j=1}\frac{1}{\sqrt{2\pi}\sigma} \exp \left(\frac{−[y_j−a \sin(\omega t_j)−b \cos(\omega t_j)]^2}{2\sigma^2} \right) $$

- where $y_i$ is the measurement (e.g., the brightness of a star) taken at time $t_i$.

- With a lof of math we do not report here, and assuming uniform priors on $a, b, \omega$, and $\sigma$ (which gives nonuniform priors on $A$ and $\phi$), the posterior distribution of parameters can be simplified to:

$$p(\omega,a,b,\sigma|{t,y}) \propto \sigma^{−N} \exp \left(\frac{−NQ}{2\sigma^2} \right)$$

%% Cell type:markdown id:e360d1c2-5506-485c-900c-4d1b021bad36 tags:

- with

$$Q= V - {2\over N} \left[ a \, I(\omega) + b \, R(\omega) - a\, b\, M(\omega) - {1 \over 2} a^2 \, S(\omega) - {1 \over 2} b^2 \,C(\omega)\right]$$

- and

$$ V = {1\over N} \sum_{j=1}^N y_j^2$$

$$       I(\omega) = \sum_{j=1}^N y_j   \sin(\omega t_j)$$

$$ R(\omega) = \sum_{j=1}^N y_j  \cos(\omega t_j)$$

$$      M(\omega) = \sum_{j=1}^N \sin(\omega t_j) \, \cos(\omega t_j)$$

$$      S(\omega) = \sum_{j=1}^N \sin^2(\omega t_j)$$

$$ C(\omega) = \sum_{j=1}^N  \cos^2(\omega t_j)$$

- *Note that I, R, M, S, C only depend on $\omega$ and the data*.

%% Cell type:markdown id:9d5ecc2c-8d4c-47bf-9828-1e993f1030ff tags:

- If $N>>1$ and we have data that extends longer than the period:

$$S(\omega) \approx C(\omega) \approx N/2 \qquad M(\omega) \ll N/2$$

- and

$$Q \approx V - {2\over N} \left[ a \, I(\omega) + b \, R(\omega)\right]  + {1 \over 2} (a^2 + b^2)$$

%% Cell type:markdown id:5ed0717e-f81a-469b-88be-04ce1e2ba63b tags:

### The posterior for many, randomly spaced, observations

***

- If we marginalize over $a$ and $b$ (as we are interested in the period):

$$  p(\omega,\sigma|\{t,y\}) \propto  \sigma^{-(N-2)} \exp \left( { - N V \over 2 \sigma^2} + { P(\omega) \over \sigma^2}       \right)$$

- with

$$P(\omega) = {1 \over N} [ I^2(\omega) + R^2(\omega)]$$

$$             V = {1\over N} \sum_{j=1}^N y_j^2$$

$$       I(\omega) = \sum_{j=1}^N y_j   \sin(\omega t_j)$$

$$ R(\omega) = \sum_{j=1}^N y_j  \cos(\omega t_j)$$

%% Cell type:markdown id:388fafb9-b03b-4a4d-a1b0-443e2ef2e8f1 tags:

-  we know the noise $\sigma$ then

$$   p(\omega|\{t,y\}, \sigma) \propto \exp \left( { P(\omega) \over \sigma^2} \right)$$

- and we now have the posterior for $\omega$!

%% Cell type:markdown id:0c1e2930-dc56-43cd-8976-bcd984345d5b tags:

## Significance of the peaks in the periodogram

***

- Let's compute the $\chi^2$ for the LS periodogram:

$$\chi^2(\omega) \equiv {1 \over \sigma^2} \sum_{j=1}^N [y_j-y(t_j)]^2 =

  {1 \over \sigma^2} \sum_{j=1}^N [y_j- a_0\, \sin(\omega t_j) - b_0 \, \cos(\omega t_j)]^2$$

- which we can simplify to:

$$\chi^2(\omega) =  \chi_0^2 \, \left[1 - {2 \over N \, V}  \, P(\omega) \right]$$

- where, again, $P(\omega)$ is the periodogram and $\chi_0^2$ is the $\chi^2$ for a model with $y(t)$=constant:

$$  \chi_0^2 = {1 \over \sigma^2} \sum_{j=1}^N y_j^2 = {N \, V \over \sigma^2}$$

%% Cell type:markdown id:d6c806d4-cdd0-4063-a418-7cdd7571c6a8 tags:

- We'll now renormalise the periodogram as:

$$P_{\rm LS}(\omega) = \frac{2}{N V} P(\omega),$$

- where $0 \le P_{\rm LS}(\omega) \le 1$.

- With this renormalization, the reduction in $\chi^2(\omega)$ for the harmonic model, relative to $\chi^2$ for the pure noise model, $\chi^2_0$ is:

$${\chi^2(\omega) \over \chi^2_0}=  1 - P_{LS}(\omega).$$

- To determine if our source is variable or not, we first compute $P_{\rm LS}(\omega)$ and then model the odds ratio for our variability model vs. a no-variability model.

- If our variability model is "correct", then the peak of $P(\omega)$ gives the best $\omega$ and the $\chi^2$ at $\omega = \omega_0$ is $N$.

%% Cell type:markdown id:c7b31e42-6ff3-4106-bc40-0304f26893cf tags:

- If the true frequency is $\omega_0$ then the maximum peak in the periodogram should have a height:

$$P(\omega_0) = {N \over 4} (a_0^2 + b_0^2)$$

- and standard deviation:

$$      \sigma_P(\omega_0)  = {\sqrt{2} \over 2} \, \sigma^2.$$

%% Cell type:markdown id:2a98318f tags:

### Credits

***

This notebook contains material obtained from https://github.com/gnarayan/ast596_2023_Spring.

%% Cell type:markdown id:05e93b1d tags:

## Course Flow

<table>

  <tr>

    <td>Previous lecture</td>

    <td>Next lecture</td>

  </tr>

  <tr>

    <td><a href="Lecture-Lomb-Scargle.ipynb">Spectral analysis</a></td>

    <td><a href="Lecture-Lomb-Scargle.ipynb">Spectral analysis</a></td>

  </tr>

 </table>

%% Cell type:markdown id:591bd355 tags:

**Copyright**

This notebook is provided as [Open Educational Resource](https://en.wikipedia.org/wiki/Open_educational_resources). Feel free to use the notebook for your own purposes. The text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/), the code of the examples, unless obtained from other properly quoted sources, under the [MIT license](https://opensource.org/licenses/MIT). Please attribute the work as follows: *Stefano Covino, Time Domain Astrophysics - Lecture notes featuring computational examples, 2025*.

%% Cell type:code id:226679b8-f75b-456d-ab17-6043f9066819 tags:

``` julia

```

This notebook is provided as [Open Educational Resource](https://en.wikipedia.org/wiki/Open_educational_resources). Feel free to use the notebook for your own purposes. The text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/), the code of the examples, unless obtained from other properly quoted sources, under the [MIT license](https://opensource.org/licenses/MIT). Please attribute the work as follows: *Stefano Covino, Time Domain Astrophysics - Lecture notes featuring computational examples, 2026*.

# This file is machine-generated - editing it directly is not advised

julia_version = "1.12.3"

julia_version = "1.12.5"

manifest_format = "2.0"

project_hash = "31cf997e22990a4cf1332dedb6a002a5d3a05816"

project_hash = "8a66b963c064ca214f5715a98f5ae2d8ea664b09"

[[deps.AbstractFFTs]]

deps = ["LinearAlgebra"]

version = "0.1.6"

[[deps.FFMPEG_jll]]

deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "PCRE2_jll", "Zlib_jll", "libaom_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]

deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "PCRE2_jll", "Zlib_jll", "libaom_jll", "libass_jll", "libfdk_aac_jll", "libva_jll", "libvorbis_jll", "x264_jll", "x265_jll"]

git-tree-sha1 = "01ba9d15e9eae375dc1eb9589df76b3572acd3f2"

uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"

version = "8.0.1+0"

[[deps.MozillaCACerts_jll]]

uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

version = "2025.5.20"

version = "2025.11.4"

[[deps.MsgPack]]

deps = ["Serialization"]

uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"

version = "1.3.7+0"

[[deps.Xorg_libXfixes_jll]]

deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libX11_jll"]

git-tree-sha1 = "75e00946e43621e09d431d9b95818ee751e6b2ef"

uuid = "d091e8ba-531a-589c-9de9-94069b037ed8"

version = "6.0.2+0"

[[deps.Xorg_libXrender_jll]]

deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libX11_jll"]

git-tree-sha1 = "7ed9347888fac59a618302ee38216dd0379c480d"

uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa"

version = "0.9.12+0"

[[deps.Xorg_libpciaccess_jll]]

deps = ["Artifacts", "JLLWrappers", "Libdl", "Zlib_jll"]

git-tree-sha1 = "4909eb8f1cbf6bd4b1c30dd18b2ead9019ef2fad"

uuid = "a65dc6b1-eb27-53a1-bb3e-dea574b5389e"

version = "0.18.1+0"

[[deps.Xorg_libxcb_jll]]

deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libXau_jll", "Xorg_libXdmcp_jll"]

git-tree-sha1 = "bfcaf7ec088eaba362093393fe11aa141fa15422"

uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

version = "5.15.0+0"

[[deps.libdrm_jll]]

deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libpciaccess_jll"]

git-tree-sha1 = "63aac0bcb0b582e11bad965cef4a689905456c03"

uuid = "8e53e030-5e6c-5a89-a30b-be5b7263a166"

version = "2.4.125+1"

[[deps.libfdk_aac_jll]]

deps = ["Artifacts", "JLLWrappers", "Libdl"]

git-tree-sha1 = "646634dd19587a56ee2f1199563ec056c5f228df"

uuid = "a9144af2-ca23-56d9-984f-0d03f7b5ccf8"

version = "1.0.21+0"

[[deps.libva_jll]]

deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libX11_jll", "Xorg_libXext_jll", "Xorg_libXfixes_jll", "libdrm_jll"]

git-tree-sha1 = "7dbf96baae3310fe2fa0df0ccbb3c6288d5816c9"

uuid = "9a156e7d-b971-5f62-b2c9-67348b8fb97c"

version = "2.23.0+0"

[[deps.libvorbis_jll]]

deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll"]

git-tree-sha1 = "11e1772e7f3cc987e9d3de991dd4f6b2602663a5"

CommonMark = "a80b9123-70ca-4bc0-993e-6e3bcb318db6"

DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"

IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a"

LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"

LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

LombScargle = "fc60dff9-86e7-5f2f-a8a0-edeadbb75bd9"

PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"