notebook.community

Edit and run 





In [91]:



    


%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn
plt.rcParams['axes.labelsize'] = 22

Estimating the underlying distributtion from correlations in $\Delta\dot{\nu}$ and $|\dot{\nu}|$

In the supplementary material of Lyne (2010) the authors noted correlation between the $\log$ values of the fractional amplitude $\Delta\dot\nu$ (which they interpreted as the switching magnitude) and the pulsars spin-down rate $|\dot{\nu}|$.

Reproduce the correlation

First let's reproduce the result, here is the data from Table 1:





In [ ]:



    


nu = [1.229, 1.616, 1.543, 10.4, 1.3, 2.579, 2.622, 1.410, 2.631, 1.672,
      3.952, 1.524, 0.983, 2.469, 3.256, 2.322, 5.996, 3.728]

nudot = [-21.25, -12.05, -11.76, -135.36, -88.31, -11.84, -7.5, -1.75, -1.18,
         -17.7, -3.59, -22.75, -5.33, -365.68, -58.85, -73.96, -604.36, -58.64]

Deltanudot_over_nudot_percent = [44.9, 13.28, 6.8, 5.91, 3.28, 2.52, 2.0, 1.71,
                         1.69, 0.85, 0.84, 0.79, 0.75, 0.71, 0.68, 0.68, 0.66,
                         0.31]

Now turn into into S.I units and calculate the $\Delta\dot{\nu}$ from the fractional size





In [3]:



    


Deltanudot_over_nudot = np.array(Deltanudot_over_nudot_percent) / 100.0

nudot = 1e-15 * np.array(nudot)
Deltanudot = Deltanudot_over_nudot * np.abs(nudot)

Plot the data





In [93]:



    


xobs = np.log(np.abs(nudot))
yobs = np.log(np.abs(Deltanudot))

p = np.corrcoef(xobs, yobs)[0, 1]
print("Correlation coefficient of {:1.3f}".format(p))

fig, ax = plt.subplots()
ax.plot(xobs, yobs, "o")
ax.set_xlabel(r"$\log_{10}(|\dot{\nu}|)$")
ax.set_ylabel(r"$\log_{10}(\Delta\dot{\nu})$")
plt.show()



    


















    






Correlation coefficient of 0.768

This seems to agree with what was found by Lyne (2010), although they had an extra 51 data points from (S2).

Precession intepretation

If the variations are caused not by switching, but instead by precession, then the time-dependent spin-down modulatoin is given by

$$ \Delta\dot{\nu}(t) = \frac{1}{2\pi\tau_A P} \left(2\theta\cot\chi\sin\psi(t) - \frac{\theta^{2}}{2}\cos2\psi(t)\right) $$

The $\Delta\dot{\nu}$ considered in the previous plot are the peak-to-peak values. Solving for the extrema we calculate (see the end of this notebook) that the precession predicts the peak-to-peak amplitude to be

$$ |\Delta\dot{\nu}| = \frac{1}{2\pi\tau_A P} \left(\left|2\theta\cot\chi\right| + \theta^{2} + \cot^{2}\chi\right) $$

Quick numerical check:

Here is a quick numerical experiment to ensure that this equation holds up





In [165]:



    


def f(t, psidot, theta, chi):
    psi = t *psidot
    return 2*theta*np.cos(chi)/np.sin(chi)*np.sin(psi) - 0.5*theta**2*np.cos(2*psi)

times = np.linspace(0, 10, 1000)
psidot = 2
theta = 0.1
chi = 1.58
y = f(times, psidot, theta, chi)
plt.plot(times, y)
plt.show()

print y.max() - y.min()
print abs(2*theta*np.cos(chi)/np.sin(chi)) + theta**2 + (np.cos(chi)/np.sin(chi))**2



    


















    

















    






0.0119254746931
0.0119254990018

Predicting the underlying $\theta$ and $\chi~$

Let's first rewrite the prediction for $\Delta\dot{\nu}$ as a function of $\nu$ as follows:

\begin{align} |\Delta\dot{\nu}| & = \frac{1}{2\pi\tau_A P} \left(\left|2\theta\cot\chi\right| + \theta^{2} + \cot^{2}\chi\right) \\ & = \frac{1}{2\pi}\frac{|\dot{\nu}|}{\nu}\nu \left(\left|2\theta\cot\chi\right| + \theta^{2} + \cot^{2}\chi\right) \\ & = \frac{1}{2\pi}|\dot{\nu}| \left(\left|2\theta\cot\chi\right| + \theta^{2} + \cot^{2}\chi\right) \\ \end{align}

Then the $\log$ of these quantities should satisfy: \begin{align} \log|\Delta\dot{\nu}| = \log|\dot{\nu}|

  • \log\left(\frac{1}{2\pi} 
          \left(\left|2\theta\cot\chi\right| + \theta^{2} + \cot^{2}\chi\right)\right)
    \end{align}

Now we can fit this model to the observed values to derive the underlying distribution of $\theta$ and $\chi$.

Fitting the distribution

To fit the distribution we use linear regression with the offset defined as above. We will use a uniform prior on $\theta$ and $\chi$ over $[0, \pi]$ and a uniform prior on the noise-parameter $\epsilon$ of $[0, 10]$.





In [195]:



    


import pymc3 as pm

with pm.Model() as model:
    theta = pm.Uniform('theta', lower=0.0, upper=np.pi)
    chi = pm.Uniform('chi', lower=0.0, upper=np.pi)
    eps = pm.Uniform('eps', lower=0.0, upper=10)
    
    mu = xobs + pm.log((pm.sqrt((2*theta*pm.cos(chi)/pm.sin(chi))**2)
                        + theta**2 + (pm.cos(chi)/pm.sin(chi))**2)/(2*np.pi))/np.log(10)
    
    # Likelihood (sampling distribution) of observations
    Y_obs = pm.Normal('Y_obs', mu=mu, sd=eps, observed=yobs)

Run the MCMC simulations





In [198]:



    


nsamples = 100000
nburn = 10000
with model:
    step = pm.Metropolis()
    trace = pm.sample(nsamples, step, progressbar=True)



    


















    






 [-----------------100%-----------------] 100000 of 100000 complete in 40.2 sec

Plot the estimated posterior distribution and traces:





In [199]:



    


pm.traceplot(trace[nburn:])
plt.show()

Corner plot





In [210]:



    


from triangle import corner
theta_chain = trace[nburn:]['theta']
chi_chain = trace[nburn:]['chi']
eps_chain = trace[nburn:]['eps']
data = np.vstack((theta_chain, chi_chain, eps_chain)).T
corner(data)
plt.show()



    


















    






(90000, 3)

Plot 100 realisations with the data





In [194]:



    


npp = 100
repeats = 100
fig, ax = plt.subplots(figsize=(10, 10))
xfit = np.linspace(xobs.min(), xobs.max(), 100)

for i  in np.arange(0, npp, 1):
    s = trace[np.random.randint(nburn, nsamples)]
    chi = s['chi']
    theta = s['theta']
    cotchi = np.cos(chi)/np.sin(chi)
    eps = s['eps']
    yfit = xfit + np.log10((np.abs(2*theta*cotchi) + theta**2 + cotchi**2)/(2*np.pi))
    ax.plot(xfit, yfit, "-", lw=0.5, color='k')
    ax.fill_between(xfit, yfit-eps, yfit+eps, alpha=0.01, color='k')

ax.plot(xobs, yobs, "o", color='r', markersize=10)
ax.set_xlabel(r"$\log_{10}(|\dot{\nu}|)$")
ax.set_ylabel(r"$\log_{10}(\Delta\dot{\nu})$")

plt.show()

Derivation of peak-to-peak amplitude

If we define

$$ y = A \sin(x) - B\cos(2x) $$

It can be shown that this has four extrema in $[0, 2\pi]$ which are

$$ x = \sin^{-1}\left(\frac{-A}{4B}\right), \; \pi - \sin^{-1}\left(\frac{-A}{4B}\right), \; \frac{\pi}{2}, \; \frac{3\pi}{2} $$

Let's check this (and help decide which is which):





In [149]:



    


def y(x, A, B):
    return A*np.sin(x) - B*np.cos(2*x)

x = np.linspace(0, 2*np.pi, 1000)
A = 1
B = 0.8

fig, ax = plt.subplots()
yplot = y(x, A, B)
ax.plot(x, yplot)

ax.set_xticks([0, np.pi/2, np.pi, 3*np.pi/2, 2*np.pi])
ax.set_xticklabels(["0", "$\pi/2$", "$\pi$", "$3\pi/2$", "$2\pi$"])

actual_dy = yplot.max() - yplot.min()
predic_dy = abs(A) + 2*B + A**2 / (8*B)
print "Actual dy = {:1.3f}, predicted df = {:1.3f}".format(actual_dy, predic_dy)
plt.show()



    


















    






Actual dy = 2.756, predicted df = 2.756

From this it is clear the last two solutions correspond to the maxima while the other two correspond to the minima. Then let's calculate the peak-to-peak amplitude. The two minima are always equal, so we have two peak-to-peak amplitudes to consider:

$$ \Delta y_1 = A + 2B + \frac{A^{2}}{8B} $$

and

$$ \Delta y_2 = -A + 2B + \frac{A^{2}}{8B} $$

The second value will be larger if $A<0$, so we can write the general peak-to-peak amplitude as

$$ \Delta y = |A| + 2B + \frac{A^{2}}{8B} $$

Provided $B > 0$

Content source: ga7g08/ga7g08.github.io

Similar notebooks:

notebook.community | gallery | about