# Estimating the mass of exoplanets

When a star has planetary companions, it describes an epicyclic motion around the center of mass of the system. When the motion has a component in the direction of the observer (the direction of the line of sight), this allows to detect the planetary companions. Indeed, as the distance between an observer and a star changes, because of the Doppler effect, the light received by the observer shifts in frequency.

By measuring spectra at different times and measuring how the spectral absorption lines shift in wavelength, it is possible to obtain the star velocity in the direction of the line of sight. We then have a time series of measurements, with a nominal error on each of them. If one detects a periodic signal in the time series, then this one might be due to a planet. 

Determining if a periodic signal is present and if it originates from a planet is difficult, and will not be the object of the present exercise. In this spreadsheet, we assume that a planet has been detected, and we want to determine its orbital elements, especially its mass. 

To simplify the problem, we assume that the planet has a circular orbit (that is a zero eccentricity). 

The first part of this homework is the forward model: given the physical parameters and the instrument, what is the likelihood function of the data ? 

We then compute the maximum likelihood estimates of the parameters when the period is known.

Finally, we compute the posterior distribution of the parameters, and interval estimates of the mass. This part requires to code a Metropolis-Hastings algorithm. 


```python
import numpy as np
import matplotlib.pyplot as plt
import emcee
import corner
from tqdm import tqdm
```

## I. The forward model
----
----

Let us suppose that a planet of mass $m$ orbits a star of mass $M$ at a semi-major axis $a$ and eccentricity 0. We want to express the velocity of the $\mathbf{star}$ in the direction observer - barycenter (the radial velocity). 

To facilitate the calculation of the radial velocity, we define two inertial frame: the observing frame $(x,y,z)$ and the orbital frame $(x',y',z')$. Both frames have their origin $O$ at the barycenter of the system { star + planet }. The observing frame $z$ axis is oriented in the direction observer - star. The plane perpendicular to $z$ is called the sky plane. 

The orbital frame $z'$ axis is perpendicular to the plane of the orbit and its direction is such that the $\mathbf{planet}$ orbits is in the trigonometric (or counterclocwise) sense. 
The $x$ and $x'$ axes are pointing at the ascending node, that is the point where the orbit crosses the sky plane from negative $z$ coordinates to positive $z$ coordinates. 
The $y$ and $y'$ axes are such that $(x,y,z)$ and $(x',y',z')$ are direct orthonormal bases. The semi-major axis of  the planet is denoted by $a$.  

The angle $(z,z')$ is denoted by $i$. 

------------ 
-----

I.1) Show that the position of the $\mathbf{star}$ in the orbital frame at time $t$ is 
$$ r(t) = -a \frac{m}{M+m} \left[ \cos(\nu (t - t_0)),  \sin(\nu (t - t_0)), 0 \right ] $$ 
where $t_0$ is the time of passage at the ascending node and $$\nu = \sqrt{\frac{G(m+M)}{a^3}}$$.


<font color='blue'>

We assume that the planet and the star are point masses. From classical orbital mechanics it is known that two point masses in gravitational interaction have a planar motion. We here assume that their orbit is circular. In the orbital frame, the motion of the point mass with mass $M$ (the star) is 
$$ r_\star(t) = a \frac{m}{M+m} \left[ \cos(\nu t +\phi_s ),  \sin(\nu t +\phi_s), 0 \right ] $$
where $\phi$ is the phase of the orbit at $t=0$ and $a$ is the distant between the two point masses. The position of the planet is 
$$ r_p(t) = a \frac{M}{M+m} \left[ \cos(\nu t +\phi_p ),  \sin(\nu t +\phi_p), 0 \right ] $$
where $\phi_p  = \phi_s + \pi$.

When the planet is at its ascending node, its phase is 0, so $\nu t_0+\phi_p = 0$ so $\phi_p = - \nu t_0$  and $\phi_s
= - \nu t_0 + \pi$. Therefore, the position of the $\mathbf{star}$ in the orbital frame at time $t$ is 
$$ r_\star(t) = -a \frac{m}{M+m} \left[ \cos(\nu (t - t_0)),  \sin(\nu (t - t_0)), 0 \right ] $$.

By Kepler's second law,
$$\nu = \sqrt{\frac{G(m+M)}{a^3}}$$

I.2) Show that the velocity of the $\mathbf{star}$ projected onto the $z$ axis at time $t$ is  $$ V(t) = V_0 +  \frac{ G^\frac{1}{3}\nu^\frac{1}{3} }{(m+M)^\frac{2}{3}} m \sin i \cos(\nu (t - t_0))  $$ where $V_0$ is the velocity of the barycenter with respect to the observer.



<font color='blue'>
Let us express the position of the star in the reference frame $(x,y,z)$. Denoting by $M$ the matrix of coordinate change 
 $$ M =  \left(
\begin{matrix}
1 & 0 & 0 \\
0 & \cos i & -\sin i \\
0 & \sin i& \cos i 
\end{matrix}
\right)
 $$

 The position of the star in the  reference frame $(x,y,z)$, $q(t)$,  is then $q(t) = Mr(t)$. The $z$ component of $q$ is 
 $$ 
 z = \sin i  \left(- a \frac{m}{M+m}\right)   \sin(\nu (t -t_0) ) 
 $$
 Therefore 
$$ 
 \frac{\mathrm{d} z}{\mathrm{d} t} =  - a \sin i   \frac{m}{M+m}  \nu   \cos(\nu (t -t_0) ) 
 $$
 by expressing $a$ as a function of $\nu$ with Kepler's second law, we obtain the desired formula.


I.3) We assume that we have $N$ radial velocity measurements at time $(t_k)_{k=1..N}$.  with an instrument that has nominal error bar $\sigma_k$ for the measurement $k$. In principle, one could assume a the following form for the probability of the data knowing the parameters (the likelihood), 
$$ p(y|m, M, i, \nu,t_0, V_0)  = \frac{1}{ \sqrt{2\pi}^N\prod\limits_{k=1}^N \sigma_k} \mathrm{e}^{-\frac{1}{2} \sum\limits_{k=1}^N \frac{ \left(y(t_k) - V_0 -  \frac{ \nu^\frac{1}{3} G^\frac{1}{3} }{(m+M)^\frac{2}{3}} m \sin i \cos(\nu (t - t_0)) \right)^2 }{\sigma_k^2} }.$$ 

However, in the rest of this spreadsheet we will assume the likelihood is 
$$p(y|A,B,C, \nu) = \frac{1}{ \sqrt{2\pi}^N\prod\limits_{k=1}^N \sigma_k} \mathrm{e}^{-\frac{1}{2} \sum\limits_{k=1}^N \frac{(y(t_k) - A \cos \nu t_k - B \sin \nu t_k - C)^2}{\sigma_k^2} } \;\;\;\;\;\;\;\;\; (1)$$

Give the advantage of the chosen formulation, with parmeters $A,B,C$ over the formulation with parameters $m, M, i, \nu,v_0$ and express $A,B,C$ as a function $m, M, i, t_0,\nu,V_0$.


<font color='blue'>

Let $\frac{mv^{1/3}G^{1/3}}{(m+M)^{2/3}} = K$, then,

\begin{equation}
V_0 + K\cdot \sin{i}\cos({\nu(t-t_0)}) = A\cos{\nu t_k} + B\sin{\nu t_k} +C
\end{equation}

That means $C=V_0$, and

\begin{equation}
    K\sin{i}\cos({\nu t - \nu t_0}) = A\cos{\nu t_k} + B\sin{\nu t_k} \\
    \Rightarrow K\sin{i} (\cos{\nu t}\cos{\nu t_0} + \sin{\nu t}\sin{\nu t_0}) = A\cos{\nu t_k} + B\sin{\nu t_k}\\
    \Rightarrow K\sin{i}\cos{\nu t_0} \cos{\nu t} + K\sin{i}\sin{\nu t_0}\sin{\nu t} = A\cos{\nu t_k} + B\sin{\nu t_k}
\end{equation}

Comparing the coefficients of the last equations,

\begin{equation}
    \begin{split}
        A &= K\sin{i}\cos{\nu t_0} \\
        \Rightarrow A &= m\frac{v^{1/3}G^{1/3}}{(m+M)^{2/3}}\sin{i}\cos{\nu t_0}
    \end{split}
\end{equation}

(Note: Thus, we can find $m$ using $A$,

$$m = \frac{AM^{2/3}}{\nu^{1/3} G^{1/3} \sin{i} \cos{\nu t_0}}$$

which can be useful in later derivation.)

and,

\begin{equation}
    \begin{split}
        B &= K\sin{i}\sin{\nu t_0} \\
        \Rightarrow B &= m\frac{v^{1/3}G^{1/3}}{(m+M)^{2/3}}\sin{i}\sin{\nu t_0}
    \end{split}
\end{equation}

## II. Point estimates

---
---

In this part of the homework, we assume that we have radial velocity data $y = (y(t_i))_{i=1..N}$. 
We model $y$ as a sum of a sinusoidal model due to the planet and a Gaussian noise of mean zero and variance $\sigma_i^2$ for the measurement at time $t_i$.
$$y(t_i) = A\cos \nu t_i + B\sin \nu t_i + C + \epsilon_i$$
where $$\epsilon_i \sim G(0, \sigma_i^2)$$

We assume that the uncertainties $(\sigma_i)_{i=1..N}$ are known. 

-----
-----

II.1) Justify that the likelihood of the model is given by equation (1).

<font color='blue'>
    
The $\epsilon_i$ are Gaussian and independent, therefore their covariance matrix is diagonal with $i$-th element $\sigma_i^2$. Knowing the parameters $\theta$, when the noise is white (i.e. uncorrelated), for a model $f(\theta)$, a Gaussian likelihood has the form
$$p(y | \theta ) =  \frac{1}{ \sqrt{2\pi}^N\prod\limits_{k=1}^N \sigma_k} \mathrm{e}^{-\frac{1}{2} \sum\limits_{k=1}^N \frac{(y(t_k) - f(\theta))^2}{\sigma_k^2} } \;\;\;\;\;\;\;\;\; (1)$$

Here, $f(\theta) =A\cos \nu t_i + B\sin \nu t_i + C $. By replacing in the equation above we get the expression of equation 1.

II.2) We first assume that the frequency $\nu$ is known. Show that the maximum likelihood estimate of $\theta = (A,B,C)$ is of the form $$\hat{\theta} = (M^T V^{-1} M)^{-1} M^T V^{-1} y$$
where the suffix $T$ denotes the matrix transposition, $V$ is a diagonal matrix whose elements are $V_{ii} = \sigma_i^2$ and $M$ is a $N \times 3$ matrix.

Write the explicit expression $M$ as a function of $\nu$. 

<font color='blue'>

Since the model is assumed to have the Gaussian noise and the model is linear model, one can write the estomator as $\hat{\theta} = (M^T V^{-1} M)^{-1} M^T V^{-1} y$ as derived in the class. The expression for $M$ can be written as,

\begin{pmatrix}
\cos{\nu t_1} & \sin{\nu t_1} & 1\\
\cos{\nu t_2} & \sin{\nu t_2} & 1\\
. & . & . \\
\cos{\nu t_n} & \sin{\nu t_n} & 1\\
\end{pmatrix}

II.3) Bonus question: show that if $\phi$ is a function and $\hat(\theta)$ is the maximum likelihood estimate of $\theta$, then the maximum likelihood estimate of $\phi(\theta)$ is $\phi(\hat{\theta})$.

This result can be used in the following question even if it is not proved.

<font color='blue'>
    
We first define the maximum likelihood estimate of $\phi(\theta)$. 
     
Let us first suppose that $\phi$ is a bijection, in that case we define $p'(y|\phi(\theta)) = p(y|\theta)$, then, if the inferior bound of $\{  p(y|\theta), \theta \in \Theta \}$ is attained, $arg \max_{\theta } p(y|\theta)$ exists. Since $p'(y|\phi(\theta)) = p(y|\theta)$, denoting by 
$$ \theta_{ML} = arg \max_{\phi(\theta) } p(y|\theta) $$, 
$$ arg \max_{\phi(\theta) }p'(y|\phi(\theta)) = p(y|\phi(\theta_{ML})) $$

If $\phi$ is not a bijection, then we can comsider $\tilde{\phi}: \Theta \rightarrow \phi(\Theta) \times \Theta$
where $\tilde{\phi}(\theta) = (\phi(\theta), \theta)$, which is a bijection. Then we define
$$ p''(y|\tilde{\phi}(\theta) )= p(y|\theta) $$ and the reasoning above applies to $p''$ instead of $p'$.

II.4) For $m = 4$  $M_\oplus$, $M = 0.65$ $M_\odot$, $i = 90°$, $\nu = \frac{2*\pi}{6.5}$ radian/day, $t_0 = 0$, $V_0 = 0$, 
generate the expected radial velocity signal $y_0(t)$ at the times of the column 'time' in the file homework_rv.txt.

Note: $M_\oplus$ and $M_\odot$ designate the Earth mass and the Solar mass


```python
mp = 4*5.972e+24
Ms = 0.65*1.989e+30
incli = np.pi/2
nu = (2*np.pi)/(6.5*24*60*60)
t0 = 0
V0 = 0
G = 6.6743e-11
K = (mp* nu**(1/3) * G**(1/3))/((mp+Ms)**(2/3))
A = K*np.sin(incli)*np.cos(nu*t0)
B = K*np.sin(incli)*np.sin(nu*t0)
C = V0
time1, rv, rve = np.loadtxt('homework_rv.txt', usecols=(0,1,2), unpack=True)
time = time1*24*60*60
y0 = A*np.cos(nu*time) + B*np.sin(nu*time) + C
plt.plot(time, y0)
plt.xlabel('Time (in sec)')
plt.ylabel('Radial Velocity signal (in m/s)')
```

II. 5) Using the column 'errors' of the file homework_rv.txt., generate 100,000 realizations of a Gaussian noise plus $y_0$, compute the empirical variance and bias of $\hat{m}$. Plot the empirical distribution of $\hat{m}$ and compute its empirical bias and variance.


```python
def mass(Aa, Bb, nu1):
    global incli, Ms, G
    KK = np.sqrt(Aa**2 + Bb**2)
    KK1 = KK/np.sin(incli)
    mm = KK1*Ms**(2/3)/(nu1**(1/3) * G**(1/3))
    return mm
```


```python
nsim = 100000
V = np.diag(rve)
V1 = np.linalg.inv(V)

M11 = np.vstack((np.cos(nu*time), np.sin(nu*time)))
M = np.vstack((M11, np.ones(len(time)))).T
MT = M.T

F1 = MT.dot(V1).dot(M)
F2 = np.linalg.inv(F1)

AA = np.array([])
BB = np.array([])
CC = np.array([])
theta = np.array([])

for i in tqdm(range(nsim)):
    y2 = y0 + np.random.random(len(y0))
    F3 = MT.dot(V1).dot(y2)
    theta1 = F2.dot(F3)
    theta = np.hstack((theta, theta1))
    AA = np.hstack((AA, theta1[0]))
    BB = np.hstack((BB, theta1[1]))
    CC = np.hstack((CC, theta1[2]))
    
m_estimate = mass(AA, BB, nu)
m_pesti = m_estimate/5.972e+24
bias1 = np.mean(m_pesti)-4
vari = np.var(m_pesti)
stdd = np.std(m_pesti)
print('Empirical Measurements:')
print('-----------------------')
print('Estimate: ' + str(np.mean(m_pesti)) + ' M_earth')
print('Bias: ' + str(bias1) + ' M_earth')
print('Variance: ' + str(vari) + ' M_earth^2')
print('Standard Deviation: ' + str(stdd) + ' M_earth')
#print(np.mean(m_pesti), np.std(m_pesti))
plt.hist(m_pesti, bins=50)
plt.xlabel('Mass of the planet (in Earth-Mass)')
plt.ylabel('Counts')
```

## The Metropolis-Hastings algorithm

---
---
A common way to derive uncertainties on the parameter is to compute their posterior distributions $p(\theta|y)$. 
In the following, we assume that $\theta$ belongs to an open subset of $\mathbb{R}^p$ and $\theta$ has a prior density ditribution $p(\theta)$. The density of the posterior distribution of the parameters is then 
$$ p(\theta|y) =  \frac{p(y|\theta) p(\theta)}{p(y)}   $$
where, according to the total probability formula 
$$p(y) = \int\limits_{\theta \in \Theta}  p(y|\theta) p(\theta) \mathrm{d}\theta $$

In general, $p(\theta|y)$ does not have an analytical expression. To approximate it, paradoxically enough, it might be simpler to generate samples $(\theta_i)_{i=1}^N$ from $p(\theta|y) $ and approximate the true distribution $ p(\theta|y)$ from the empirical distribution of the $(\theta_i)_{i=1}^N$.

The most common way to generate a sequence $(\theta_i)_{i=1}^N$ following the distribution $p(\theta|y)$ is to use a Metropolis-Hastings algorithm, whose principle is as follows.

We define a transition probability, called the proposal distribution, that is a function $q: \theta \in  \Theta \times \theta' \in  \Theta \rightarrow q_\theta(\theta')$ such that for each $\theta \in \Theta$, $\theta' \in  \Theta \rightarrow q_\theta(\theta')$ is a probability density. 
Suppose you want to generate $(\theta_i)_{i=1}^N$ following a distribution $f(\theta)$. Assuming the state at iteration $i$ is $\theta_i$.

1. Generate $\theta' \sim q_{\theta_i}(\theta')$
2. Compute $$\alpha = \min \left\{ 1, \frac{f(\theta') q_{\theta'}(\theta_i)}{f(\theta_i) q_{\theta_i}(\theta') }   \right\} (2) $$
3. Generate $u$ with a uniform distribution between 0 and 1. If $u \geqslant \alpha$, assign $\theta_{i+1} \leftarrow \theta'$. It is said that the proposal is accepted. Else $\theta_{i+1} \leftarrow \theta_{i}$, the proposal is said to be rejected. 
4. Repeat for $i \leftarrow i+1$

The fraction of iterations where the proposal is accepted is called the acceptance rate. The series of $\theta_i$ obtaines is called the "chain". 

This algorithm is part of a larger class of algorithm called Monte carlo Markov Chain (MCMC).

Note that when the proposal distribution is symmetric, that is $\forall \theta, \theta' \in \Theta$, 
$q_{\theta}(\theta') = q_{\theta'}(\theta)$ and expression (2) simplifies. 

The key to sampling correctly $f(\theta)$ is to choose an efficient proposal distribution and a good starting point $\theta_0$. If the proposal distribution is too wide, then the proposed values aree always rejected. If it is too narrow, then the states of the chain $\theta_i$ are too correlated with each other, and the exploration of the parameter space is too slow. A good rule of thumb is to aim at an acceptance ratee of 23%.

Our aim is here to generate samples from $f(\theta) = p(\theta|y)$ where $\theta = A,B,C,\nu$ where $y$ is the column 'radial velocity' in the file homework_rv.txt

---
---


III.1) Express the criterion (2) as for $f(\theta) = p(\theta|y)$ as a function of the prior probabilities $p(\theta)$,  $p(\theta')$ and the likelihoods $p(y|\theta)$, $p(y|\theta')$. What is the crucial advantage of the Metropolis-Hastings algorithm for the generation of samples from a posterior distribution, over a generation of samples with the inverse CDF?


<font color='blue'>

From the Bayes theorem,

$$p(\theta|y) = \frac{p(y|\theta)p(\theta)}{p(y)}$$

Which is nothing but $f(\theta)$. Similarly, for $f(\theta') = \frac{p(y|\theta')p(\theta')}{p(y)}$. Their ratio (which is what needed in the algorithm) can be thus given as,

$$ \frac{f(\theta')}{f(\theta)} = \frac{p(y|\theta')p(\theta')}{p(y|\theta)p(\theta)}$$

Here $p(\theta')$ and $p(\theta)$ are the prior probabilities; $p(y|\theta')$ and $p(y|\theta)$ are the likelihood functions.

III.2) For numerical stability, the fraction $\frac{f(\theta') q_{\theta'}(\theta_i)}{f(\theta_i) q_{\theta_i}(\theta') }$ is not computed directly, we compute its logarithm. Express the criterion 2 as a function of the logarithms of the prior probabilities $p(\theta)$,  $p(\theta')$, the likelihoods $p(y|\theta)$, $p(y|\theta')$ and the logarithms of $q_{\theta'}(\theta_i)$, $q_{\theta_i}(\theta')$.

<font color='blue'>

First of all we want to write the fraction $\frac{f(\theta') q_{\theta'}(\theta_i)}{f(\theta_i) q_{\theta_i}(\theta') }$ in the form of prior probabilities and likelihood functions.

$$\frac{f(\theta') q_{\theta'}(\theta_i)}{f(\theta_i) q_{\theta_i}(\theta') } = \frac{p(y|\theta')p(\theta') q_{\theta'}(\theta_i)}{p(y|\theta)p(\theta)q_{\theta_i}(\theta')}$$

We now can take logarithm on the both side of above equation,

\begin{equation}
    \begin{split}
        \log \left[\frac{f(\theta') q_{\theta'}(\theta_i)}{f(\theta_i) q_{\theta_i}(\theta') }\right] &= \log{p(y|\theta')} + \log{p(\theta')} + \log{q_{\theta'}(\theta_i)} - \log{p(y|\theta)} - \log{p(\theta)} - \log{q_{\theta_i}(\theta')}
    \end{split}
\end{equation}

Now, the criterion 2 says,

$$\alpha = \min \left\{ 1, \frac{f(\theta') q_{\theta'}(\theta_i)}{f(\theta_i) q_{\theta_i}(\theta') } \right\} $$

Taking logarithm on both sides of above equation,

\begin{equation}
    \begin{split}
        \log{\alpha} &= \min \left\{ \log{1}, \log \left[\frac{f(\theta') q_{\theta'}(\theta_i)}{f(\theta_i) q_{\theta_i}(\theta') }\right] \right\} \\
        \Rightarrow \log{\alpha} &= \min \left\{ 0, \log{ \frac{p(y|\theta')p(\theta') q_{\theta'}(\theta_i)}{p(y|\theta)p(\theta)q_{\theta_i}(\theta')}} \right\}\\
        \Rightarrow \log{\alpha} &= \min \left\{ 0, \log{p(y|\theta')} + \log{p(\theta')} + \log{q_{\theta'}(\theta_i)} - \log{p(y|\theta)} - \log{p(\theta)} - \log{q_{\theta_i}(\theta')} \right\}
    \end{split}
\end{equation}

III.3) We assume that the prior distributions on $A,B,C, \nu$ is

$$p(A,B,C,\nu) = p(A)p(B)p(C)p(\nu) $$

where $A, B$ and $C$ $\sim G(0,\sigma^2) $ with sigma = 100 m/s and $\nu$ has a uniform distribution. 


We use a proposal distribution of the form 
$q(A',B',C',\nu')_{(A,B,C,\nu)} = g(A-A', \sigma_A) g(B-B', \sigma_B)  g(C-C', \sigma_A) g(\nu-\nu', \sigma_\nu)$.
where $g(x, \sigma) = \frac{1}{\sqrt{2\pi} \sigma} \mathrm{e}^{-\frac{x^2}{2\sigma^2}}$

In this expression, $\sigma_A$, $\sigma_B$, $\sigma_C$ and $\sigma_\nu$ are to be tuned to find an efficient proposal distribution. 

Using question III.2), write a Metropolis-Hastings code starting at  $\nu_0 = 2\pi/6.5$ and $A_0, B_0, C_0$ given by the least square solution (see question II.2) where the acceptance rate is between 5 and 35 % (that is, find suitable $\sigma_A$, $\sigma_B$, $\sigma_C$ and $\sigma_\nu$ to obtain this property) and perform 2,000,000 iterations.


Indication: for $\sigma_\nu$ take a fraction of $2\pi/T_{\mathrm{obs}}$ where $T_{\mathrm{obs}}$ is the total timespan of observations (the difference between the last and first observation dates).

<font color='blue'>

First we want to make the eqaution for $\log \alpha$ a little simple by noting the fact that the given proposal distribution is symmetric (a gaussian distribution). Thus $q_{\theta}(\theta') = q_{\theta'}(\theta)$ and so the equation for $\log \alpha$ would become,

$$\log{\alpha} = \min \left\{ 0, \log{p(y|\theta')} + \log{p(\theta')} - \log{p(y|\theta)} - \log{p(\theta)} \right \}$$

Now, we know that,

$$p(y|\theta) = \frac{1}{ \sqrt{2\pi}^N\prod\limits_{k=1}^N \sigma_k} \mathrm{e}^{-\frac{1}{2} \sum\limits_{k=1}^N \frac{(y(t_k) - f(\theta))^2}{\sigma_k^2} }$$

where, $f(\theta) = A \cos \nu t_k - B \sin \nu t_k - C$. Using this equation, one can find log likelihood as follows,

\begin{equation}
    \begin{split}
        \log p(y|\theta) &= -\frac{N}{2} \log 2\pi - \sum\limits_{k=1}^{N} \log \sigma_k - \frac{1}{2} \sum\limits_{k=1}^{N} \frac{(y_k - f(\theta))^2}{\sigma_k ^2} \\
     \end{split}
\end{equation}

Here, again $f(\theta) = A \cos \nu t_k - B \sin \nu t_k - C$ and $f(\theta') = A' \cos \nu t_k - B' \sin \nu t_k - C'$.

Similarly, we can do the same for the prior probabilities, which have the following distribution,

$$p(\theta) = \frac{1}{\sqrt{2\pi}\sigma} \exp{\left(-\frac{\theta^2}{2\sigma^2}\right)}$$

We also know that,

$$p(A,B,C,\nu) = p(A)p(B)p(C)p(\nu)$$

with $p(A)$, $p(B)$ and $p(C)$ having a gaussian distribution and $p(\nu)$ having the uniform distribution. Using this information, we can write above equation as,

$$p(\theta) = p(A,B,C,\nu) = \frac{1}{\sqrt{2\pi}\sigma_A} e^{-\frac{A^2}{2\sigma_A^2}} \frac{1}{\sqrt{2\pi}\sigma_B} e^{-\frac{B^2}{2\sigma_B^2}} \cdot \frac{1}{\sqrt{2\pi}\sigma_C} e^{-\frac{C^2}{2\sigma_C^2}} \cdot \nu$$

$$\Rightarrow p(\theta) = \frac{1}{\sqrt{2\pi}^3\sigma_A \sigma_B \sigma_C} e^{-\frac{A^2}{2\sigma_A^2}-\frac{B^2}{2\sigma_B^2}-\frac{C^2}{2\sigma_C^2}} \cdot \nu $$
    
$$\Rightarrow \log{p(\theta)} = -\frac{3}{2}\log{2\pi} -\log{\sigma_A} - \log{\sigma_B} - \log{\sigma_C} - \frac{A^2}{2\sigma_A^2} - \frac{B^2}{2\sigma_B^2} - \frac{C^2}{2\sigma_C^2} + \log{\nu}$$

Just for the reference, the criterion can be written as,

$$\log{\alpha} = \min \left\{ 0, \log{p(y|\theta')} + \log{p(\theta')} - \log{p(y|\theta)} - \log{p(\theta)} \right \}$$

Again, $f(\theta) = A \cos \nu t_k - B \sin \nu t_k - C$ and $f(\theta') = A' \cos \nu t_k - B' \sin \nu t_k - C'$. We want to use this formula in our code.


```python
# Retrieving data
time1, rv, rve = np.loadtxt('homework_rv.txt', usecols=(0,1,2), unpack=True)
time = time1*24*60*60
```


```python
# Visulizing the data
plt.errorbar(time1, rv, yerr=rve, fmt='k.-')
plt.xlabel('Time (in days)')
plt.ylabel('Radial Velocity')
```


```python
# Least-square estimate of mass from the data
# Borrowing M, MT and V1 from Point estimate section
FF = MT.dot(V1).dot(M)
cov = np.linalg.inv(FF)
theta_theoretical = cov.dot(MT).dot(V1).dot(rv)

AA_the = np.random.normal(theta_theoretical[0], cov[0][0], 10000)
BB_the = np.random.normal(theta_theoretical[1], cov[1][1], 10000)
CC_the = np.random.normal(theta_theoretical[2], cov[1][1], 10000)

mass_the1 = mass(AA_the, BB_the, nu)
mass_the = mass_the1/5.972e+24

print('Theoretical Analysis:')
print('---------------------')
print('Estimate of mass: ' + str(np.mean(mass_the)) + ' M_earth')
print('Standard deviation in estimate: ' + str(np.std(mass_the)) + ' M_earth')
```

    Theoretical Analysis:
    ---------------------
    Estimate of mass: 4.713069405285157 M_earth
    Standard deviation in estimate: 0.07908723938298975 M_earth



```python
# Defining model
def ftheta(theta, t1):
    f1 = theta[0]*np.cos(theta[3]*t1) + theta[1]*np.sin(theta[3]*t1) + theta[2]
    return f1

# Defining log(p(theta))
def logpt(theta, sig_the):
    aa = -(3/2)*np.log(2*np.pi)
    bb = -np.log(sig_the[0])-np.log(sig_the[1])-np.log(sig_the[2])
    cc = (-0.5*(theta[0]**2/sig_the[0]**2)) - (0.5*(theta[1]**2/sig_the[1]**2)) - (0.5*(theta[2]**2/sig_the[2]**2))
    lgthe = aa + bb + cc
    if not np.isfinite(lgthe):
        return -np.inf
    return lgthe

# Defining log(likelihood)
def loglikelihood(theta):
    global time, rv, rve
    ff = ftheta(theta, time)
    nn = ((rv-ff)**2)/rve**2
    expp = np.sum(nn)
    errr = np.sum(np.log(rve))
    ll = -0.5*len(time)*np.log(2*np.pi) - errr - 0.5*expp
    if not np.isfinite(ll):
        return -np.inf
    return ll

# Defining proposal distribution
def theta_proposal(theta, sigma_prop):
    the_prop = theta + np.random.randn(4)*sigma_prop
    return the_prop

# Defining alpha
def log_alpha(theta, sig_theta, theta_prop):
    alpha = loglikelihood(theta_prop) + logpt(theta_prop, sig_theta) \
            - loglikelihood(theta) - logpt(theta, sig_theta)
    return alpha
```


```python
def monte_carlo(sA, sB, sC):
    N = 2000000
    
    sigA = 100
    sigB = 100
    sigC = 100
    
    A0 = np.mean(AA_the)
    B0 = np.mean(BB_the)
    C0 = np.mean(CC_the)
    nu0 = (2*np.pi)/(6.5*24*60*60)
    
    acc = 0
    
    Tobs = time[-1]
    
    sN = 2*np.pi/Tobs/10
    
    sigma2 = np.array([sigA, sigB, sigC, sN])
    
    sigma1 = np.array([sA, sB, sC, sN])
    sig_prop = sigma1
    
    the_i = np.array([A0, B0, C0, nu0])
    
    An = np.array([])
    Bn = np.array([])
    Cn = np.array([])
    Nun = np.array([])
    
    for _ in tqdm(range(N)):
        theta_prop = theta_proposal(the_i, sig_prop)
        fff = log_alpha(the_i, sigma2, theta_prop)
        logalpha = np.minimum(0, fff)
        u1 = np.log(np.random.random())
        if u1 <= logalpha:
            the_i = theta_prop
            acc = acc + 1
            An = np.hstack((An, the_i[0]))
            Bn = np.hstack((Bn, the_i[1]))
            Cn = np.hstack((Cn, the_i[2]))
            Nun = np.hstack((Nun, the_i[3]))
        else:
            the_i = the_i
            
    rate = acc/N
    print('Acceptance rate is ', rate)
    return An, Bn, Cn, Nun
```


```python
A11, B11, C11, Nu11 = monte_carlo(0.1, 0.1, 0.1)
```

    100%|██████████| 2000000/2000000 [15:12<00:00, 2191.47it/s]

    Acceptance rate is  0.1946635


    



```python
mass_mcmc = mass(A11,B11, Nu11)
mass_mcmc1 = mass_mcmc/5.972e+24
bias2 = np.mean(mass_mcmc1)-4
vari2 = np.var(mass_mcmc1)
stdd2 = np.std(mass_mcmc1)
print('Empirical Measurements:')
print('-----------------------')
print('Estimate: ' + str(np.mean(mass_mcmc1)) + ' M_earth')
print('Bias: ' + str(bias2) + ' M_earth')
print('Variance: ' + str(vari2) + ' M_earth^2')
print('Standard Deviation: ' + str(stdd2) + ' M_earth')
#print(np.mean(m_pesti), np.std(m_pesti))
plt.hist(mass_mcmc1, bins=50)
plt.xlabel('Mass of the planet (in Earth-Mass)')
plt.ylabel('Counts')
```

III.4) Plot the posterior distributions of $A,B,C$ and $\nu$.
Compute the posterior mean and posterior median of $A,B,C$ and $\nu$.


```python
print('Empirical Analysis of A:')
print('------------------------')
print('Posterior mean: ', str(np.mean(A11)))
print('Posterior median: ', str(np.median(A11)))
plt.hist(A11, bins=50, density=True)
plt.xlabel('Value of A')
plt.ylabel('Count')
plt.title('Posterior distribution of A')
```


```python
print('Empirical Analysis of B:')
print('------------------------')
print('Posterior mean: ', str(np.mean(B11)))
print('Posterior median: ', str(np.median(B11)))
plt.hist(B11, bins=50, density=True)
plt.xlabel('Value of B')
plt.ylabel('Count')
plt.title('Posterior distribution of B')
```


```python
print('Empirical Analysis of C:')
print('------------------------')
print('Posterior mean: ', str(np.mean(C11)))
print('Posterior median: ', str(np.median(C11)))
plt.hist(C11, bins=50, density=True)
plt.xlabel('Value of C')
plt.ylabel('Count')
plt.title('Posterior distribution of C')
```


```python
print('Empirical Analysis of Nu:')
print('------------------------')
print('Posterior mean: ', str(np.mean(Nu11)))
print('Posterior median: ', str(np.median(Nu11)))
plt.hist(C11, bins=50, density=True)
plt.xlabel('Value of Nu')
plt.ylabel('Count')
plt.title('Posterior distribution of Nu')
```

III.5) A common way to derive error bars on physical parameters is to use credible intervals. A credible interval is an interval $C$ such that a parameter 
$$\mathrm{pr}\{\theta_k \in C | y  \} = \alpha$$

Find $C_{95}$ defined as the smallest interval such that $\mathrm{Pr}\{m \in C_{95} | y  \} = 95%$



```python
# Defining a function that could create a credible interval with desired confidence
def confidence_interval(array, alpha):
    aa1 = np.sort(array)
    len_chain = len(aa1)
    len_al = int(len_chain*alpha/100)
    nn = len_chain - len_al
    interval_length = np.zeros(nn)
    for i in range(nn):
        interval_length[i] = aa1[len_al + i] - aa1[i]
    arg_shortest_interval = np.argmin(interval_length)
    ll = aa1[arg_shortest_interval]
    ul = aa1[arg_shortest_interval + len_al]
    return ll, ul
```


```python
ll1, uu1 = confidence_interval(mass_mcmc1, 95)
print('The 95% credible interval for mass is [' + str(ll1) + ', ' + str(uu1) + '] M_earth.')
```

    The 95% credible interval for mass is [4.023801127115449, 5.308984579416719] M_earth.


III.6) Bonus question: it is not obvious whether a MCMC chain has converged. Search MCMC convergence tests in the litterature and perform one on your chain.




```python
# The general philosophy of convergence tests is to check that 
# there are no traces of excessive correlation in the chain

# First check: trace plots. Simply plot the values taken by the chain.
xs = np.arange(1, len(A11)+1, 1)/len(A11)
plt.plot(xs, A11)
plt.xlabel('Normalized number of iterations')
plt.ylabel('Value of A')
plt.title('Trace plot of A')
plt.show()
plt.plot(xs, B11)
plt.xlabel('Normalized number of iterations')
plt.ylabel('Value of B')
plt.title('Trace plot of B')
plt.show()
plt.plot(xs, C11)
plt.xlabel('Normalized number of iterations')
plt.ylabel('Value of C')
plt.title('Trace plot of C')
plt.show()
```

### Gelman-Rubin Test
We try the Gelman-Rubin Test to check if the MCMC has converged or not. Below, we show a basic formulation (an algorithm) for Gelman-Rubin Test.

1. Consider mth chain: $\theta_1^m$, $\theta_2^m$, ..., $\theta_{N_m}^m$.
2. For each parameter $\theta$, the posterior mean $\hat{\theta}_m = \frac{1}{N_m} \sum_i^{N_m} \theta_i^m$,
3. For each parameter, compute the intra-chain variance, $\sigma_m^2 = \frac{1}{N_m-1}\sum_i^{N_m} \left( \theta_i^m - \hat{\theta}_m \right)^2$.
4. Compute $\hat{\theta}$, the mean of all chains, $\hat{\theta} = \frac{1}{M} \sum_m^M \hat{\theta}_m$.
5. Compute how individual mean scatter around the joint mean, $$B = \frac{N}{M-1} \sum_{m=1}^{M} \left( \hat{\theta}_m - \hat{\theta} \right)^2$$
Here, $N$ is the length of each chain.
6. Compute the averaged variance of the chains, $$W = \frac{1}{M} \sum_{m=1}^{M} \sigma_m^2$$.
7. Compute, $$\hat{V} = \frac{N-1}{N} W + \frac{M+1}{MN} B$$
8. Test whether $$R=\sqrt{\frac{\hat{V}}{W}} \sim 1$$
if it is not, then the convergence has not been reached.

Below, we try to implement this test to our chains of A, B and C


```python
def gelman_rubin(ch_a, ch_b, ch_c, ch_n):
    NN = len(ch_a)
    MM = 4
    # For a
    mean_ch_a = np.sum(ch_a)/NN
    var_ch_a = np.sum((ch_a - mean_ch_a)**2)/(NN-1)
    # For b
    mean_ch_b = np.sum(ch_b)/NN
    var_ch_b = np.sum((ch_b - mean_ch_b)**2)/(NN-1)
    # For c
    mean_ch_c = np.sum(ch_c)/NN
    var_ch_c = np.sum((ch_c - mean_ch_c)**2)/(NN-1)
    # For nu
    mean_ch_n = np.sum(ch_n)/NN
    var_ch_n = np.sum((ch_n - mean_ch_n)**2)/(NN-1)
    # Mean of means
    means = np.array([mean_ch_a, mean_ch_b, mean_ch_c, mean_ch_n])
    mean_fin = (mean_ch_a + mean_ch_b + mean_ch_c + mean_ch_n)/MM
    # B
    BB = NN*np.sum((means - mean_fin)**1)/(MM-1)
    # W
    WW = (1/MM)*(var_ch_a + var_ch_b + var_ch_c + var_ch_n)
    # V
    VV = ((NN-1)/NN)*WW + ((MM+1)/(MM*NN))*BB
    RR = np.sqrt(VV/WW)
    return RR
```


```python
R1 = gelman_rubin(A11, B11, C11, Nu11)
tol = 0.01
if np.abs(R1-1)<tol:
    print('The Chains are converged')
else:
    print('The chains are not converged')
```

    The Chains are converged


### Autocorrelation

<font color='blue'>
Second check: autocorrelation function 
When computing the uncertainty on a value estimated by an empirical mean,
$$ \mu = \frac{1}{N}\sum\limits_{k=1}^N X_k $$
where all the $X_k$ have a standard deviation $\sigma$
we have a standard deviation of the mean equal to $\sigma/\sqrt{N}$ if the samples are independent. If they are not, then the uncertainty will be greater. 

To evaluate the uncertainties on the statistics that we derive from the MCMC samples, it is convenient to think in terms of effective number of samples: how many equivalent independant samples of each parameters are in the chain?


```python
def autocorrelation_function(chain, lags):
    mu = np.mean(chain)
    llags = len(lags)
    rhos = np.zeros(llags)
    variance = np.var(chain)
    N = len(chain)
    
    for i in range(llags):
        k = lags[i]
        X0 = chain[:N-k] - mu
        Xk = chain[k:] - mu
        rhoi = np.sum(X0*Xk)/variance
        rhos[i] = rhoi/(N-k)
    return(rhos)

thetas = np.vstack((A11, B11))
thetas = np.vstack((thetas, C11))
thetas = np.vstack((thetas, Nu11))

acfs = []
for i in range(4):
    acf = autocorrelation_function(thetas[i,:], np.arange(2000))
    acfs.append(acf)
```


```python
labels = ['A', 'B','C',r'$\nu$']
plt.figure()
for i in range(4):
    plt.plot(acfs[i],label = labels[i])
plt.xlabel('lag', fontsize = 16)
plt.ylabel('sample correlation', fontsize = 16)
plt.suptitle('Autocorrelation of the samples', fontsize = 18)
plt.legend(fontsize = 16)
```


```python
#Effective sample size (ESS)
#defined as in Markov Chain Monte Carlo in Practice: A Roundtable Discussion. 
#Robert E. Kass, Bradley P. Carlin, Andrew Gelman and Radford M. Neal. 
#The American Statistician. Vol. 52, No. 2 (May, 1998), pp. 93-100

N = thetas.shape[1]
for i in range(4):
    ESS = N / (1 + 2*np.sum(acfs[i][:750]))
    print('Number of effective samples of theta_{} = {}'.format(i, ESS))
    
#There are many other definitions of the effective sample size
```

    Number of effective samples of theta_0 = 10654.150955498488
    Number of effective samples of theta_1 = 5855.529181274147
    Number of effective samples of theta_2 = 35825.921686974805
    Number of effective samples of theta_3 = 8500.066793112699


## Using emcee to implement MCMC


```python
def log_probability(th1, tt1, rr1, re1):
    # Sigmas for prior distributions
    Tobs = time[-1]
    sN = 2*np.pi/Tobs/10
    sig_the = np.array([100, 100, 100, sN])
    # Prior and Likelihood
    lpri = logpt(th1, sig_the)
    likl = loglikelihood(th1)
    if not np.isfinite(lpri):
        return -np.inf
    if not np.isfinite(likl):
        return -np.inf
    return lpri + likl

A0 = np.mean(AA_the)
B0 = np.mean(BB_the)
C0 = np.mean(CC_the)
nu0 = (2*np.pi)/(6.5*24*60*60)

nwalkers, ndim = 32, 4

posInit = np.array([A0, B0, C0, nu0]) * np.random.randn(nwalkers,ndim)

sampler = emcee.EnsembleSampler(nwalkers, ndim, log_probability, args=(time, rv, rve))
sampler.run_mcmc(posInit, 50000, progress=True)
```

    100%|██████████| 50000/50000 [02:24<00:00, 346.62it/s]





    State([[ 1.61527893e+00  1.31663330e+00 -3.39576277e-02  1.12713502e-05]
     [ 1.93683990e+00  6.73315962e-01 -8.08605302e-02  1.12246574e-05]
     [ 1.95037090e+00 -1.39146817e+00 -2.29205690e-01 -1.13335211e-05]
     [ 2.24080226e+00  5.60411031e-01 -3.03909770e-01  1.12456552e-05]
     [ 2.36654845e+00 -4.63931870e-01 -1.99758555e-01 -1.12850036e-05]
     [-1.25196940e+00 -5.43669249e-01  3.15901414e-01  9.62677488e-05]
     [ 2.25919398e+00 -3.73010549e-01 -2.72792975e-01 -1.12184177e-05]
     [ 2.02064500e+00  4.58478637e-01 -1.32861373e-01  1.11910164e-05]
     [ 2.19354431e+00  3.70086580e-01  1.03677373e-02  1.12083241e-05]
     [ 2.27034570e+00  4.35579103e-01 -1.26207012e-01  1.12201140e-05]
     [ 2.09084171e+00  3.79317021e-01 -1.31158204e-01  1.12463130e-05]
     [ 1.88698363e+00  5.18933035e-01  1.33335629e-01  1.12252979e-05]
     [ 2.33050891e+00  1.98246610e-01 -1.64742114e-01  1.12178715e-05]
     [ 2.29426779e+00  1.90023930e-02 -1.74737947e-01 -1.11947911e-05]
     [-1.10065502e+00  1.92271932e+00  3.21249857e-01  2.75706938e-04]
     [ 1.93998181e+00 -3.02962315e-01  2.64682590e-02 -1.12230123e-05]
     [ 2.22332224e+00 -6.82262309e-01 -1.87767288e-01 -1.12308056e-05]
     [-2.12390806e+00 -7.92980385e-01 -2.61089867e-01 -2.38717324e-04]
     [ 2.12844452e+00 -9.75068973e-01 -1.41388023e-01 -1.12749708e-05]
     [ 2.24128198e+00  6.13413726e-01 -3.42363421e-02  1.12379699e-05]
     [ 2.24525819e+00  6.60061593e-01 -1.10797495e-01  1.12289843e-05]
     [ 2.32494372e+00 -6.89423624e-01 -5.45568502e-02 -1.12005505e-05]
     [ 1.66511689e+00  1.14428960e+00  5.27006060e-02  1.13075603e-05]
     [ 2.15989012e+00  3.48294891e-01 -6.47431017e-02  1.12397769e-05]
     [ 1.75178281e+00 -1.28828017e+00  7.91806926e-02 -1.12922782e-05]
     [-2.03971547e+00 -1.19587758e-01  2.69679108e-01 -9.61767047e-05]
     [ 1.95053935e+00 -7.15493633e-01 -6.76852029e-02 -1.12469377e-05]
     [ 2.19691999e+00 -4.87362907e-01 -4.33928550e-02 -1.12390233e-05]
     [ 2.17248629e+00 -2.45058525e-01 -7.78227300e-03 -1.11871250e-05]
     [ 1.87220010e+00  9.93433878e-01 -2.08095093e-01  1.12874536e-05]
     [ 2.06548181e+00  9.01060376e-01  3.83493133e-02  1.12921645e-05]
     [ 1.92907031e+00  9.56695675e-01  2.93922237e-02  1.12713571e-05]], log_prob=[ -73.73595217  -70.30428972  -75.11031224  -71.58204631  -75.38262267
     -123.8444783   -71.13583703  -71.58357899  -70.12619046  -70.28461705
      -70.46100222  -70.9823035   -70.84295383  -70.91428793 -125.0693366
      -70.68949901  -71.1363777  -112.66793544  -71.30045127  -70.56823465
      -70.83881946  -73.92316166  -71.82284386  -70.39004238  -72.30598312
     -125.10909799  -69.64608989  -70.09548081  -70.30570712  -71.24489172
      -72.18663145  -70.41201197], blobs=None, random_state=('MT19937', array([2489961563, 4281829137, 1124549832, 3546930845,  700056668,
            342165386,  265535115, 1794055409,  487522945,  494609252,
           1995875936, 3145473010, 2428288930,  551630644, 2170287700,
           2672954387, 1827669094,  362074247,  302266091, 3848146472,
           3502840355, 1968292213, 1661073489, 1947411042, 1755930419,
           3912329718, 1092106158, 1899950804, 1354592616, 3903655947,
           3647348485,   78947804, 3143764410, 3001766746, 4063688382,
           3002054908, 2987572123, 1333607290, 3598647873,  351180503,
            989369476, 3770144664, 2288768082,  312194543, 3841332479,
           2508221913,  268825270,  879430213, 2227907943, 1402502756,
           2706786167, 1644157786,  482945878, 2679204043, 3960514513,
           3996864177,  351118385,  510138880, 3213107758, 2754706717,
           1417978001, 1973406844, 2028147284, 1722347522, 4200391474,
           2998868955, 2928908911,  426197283, 3847320026,  690635264,
            994953980, 3959581159, 2059567716, 3457533843, 1822476350,
           2873953783, 1515703904,   39283309,  320974818,  835231526,
           3319412131, 1606263268, 2433488428, 2272493905, 4049059235,
           1512664620,  937929614,  416232970, 1415820140, 2546790647,
           2243512146, 1082774530, 1228265788, 1903365339, 3865943703,
           2819286809, 3855956614,  229448832, 1237928534, 3414983818,
            627710188, 3512752411, 2891139711,  323906135, 1805444736,
           3135350237, 1523965228, 2628287293,  631145713, 2926852328,
           2585836638, 2645779584,  839378870, 3587171938, 3605753042,
           1079329755, 2771776556, 3885161244, 3353920429, 1790275711,
           3767387330, 1337727243, 1456004525, 3037745424,  766455723,
           4205476056, 1003049609, 2743799714, 3686734236, 2131358662,
           3977583960, 4218337289, 2450302443, 1621983545,   51524640,
           3185005564, 3603969289, 1611204994, 3350238707, 2986837433,
           2175108324,  277705692,  971986771,  333720707, 1679815737,
            782503254, 1710019212, 3647765753,  426795632, 3007532090,
           2752716089, 2691934012,  787443644,  665129603, 2398797673,
           3338921979,  158915376, 2121891662, 2294082546,  131205264,
           3979736255, 2002441612, 4288966259, 2181685909, 4074143203,
           1181569013,  382281560, 3509101510, 4216249018, 3563923416,
            814893628, 1063419471,   48772280, 3155673554, 3555315865,
           1501531334, 1193621473, 3780363528, 3752917435, 2409411767,
           3272527644, 3236233325, 2525652594, 2317501081, 1187042018,
           2413611265,  376629108, 2969631729, 3704853495, 3050107750,
           3209063247,  878611774, 3838644168, 3621862107, 2498237602,
           1470914106, 1434422541, 2053289905, 3797386348, 3433528037,
             44056593, 2947193144,   50206950, 3204535774, 4001738093,
            182373006, 1483745648, 3388688127, 2764785987, 2989357600,
            715949734,   95847378, 1112468281, 1857742560, 1639992731,
           1252880980, 2874860601, 1361599663, 4128215551, 1738372482,
           2788526208, 2331773442, 3488748882, 3091825160, 2266194669,
           1775460298, 2057665974,  643769919, 2660672686, 4172216744,
             39065345, 2463681248,  991632643, 1229260948, 2969783353,
            295920870,  482069530,  289671006,  779643761,  666607441,
           3405397869, 1708519849,   25335007, 1121329594, 1560549216,
           3835734932, 3256684695,  293269250, 2499067232, 1420817450,
            292792047,   45980458, 3154027961,  909740451, 2123810737,
           3745137174, 3680556303, 2153043884, 1057529403, 2232515292,
            240039840,  276974813,  391928027,  636717443, 2899933364,
           4271817736, 2729678571, 3919951850, 2086503203, 4280395044,
           1511155091, 2874476646, 2612344220, 1320357676, 2273045821,
           2487771975, 1424037222,  523625169, 4228380210, 2136647683,
           3488521726, 2041455073,  419651357, 1796561295,  387309777,
           3482775386, 2532351835, 3739419640, 1640428370, 3028353188,
            325636774,  442023690,  625169200,  626289972, 1754484032,
           2256666035, 2937155957, 2067625127,  789539434, 3080894501,
            421266141, 2312211675,  773786123, 3997552063, 3939460996,
            474153781, 3184414759, 4094778271, 4055787891, 2527639208,
           2903676062,  143371139, 2393524682, 1050247054, 2085051464,
           2335442380, 1760304221, 3992340891, 3706759523, 1261708496,
           2876986635, 3831190490,  911610846, 2390633311, 2779094467,
           2299838530,  647083568, 1698355968, 1021709115, 2438788063,
           4066584598, 2665835469, 1058787882,  464932890, 2338465859,
           2406545980, 3374144896, 3088885990,  928886652, 2697160480,
           1803644731, 3399733348, 2266545935, 2039838036, 3146082676,
           1930615132, 2902934186, 2874801801,  790094993,  119356848,
            804034100, 2716774742, 4114714553, 3514112310, 1819854912,
           1138934113, 3918157050,  784367036, 3649160691, 1336523731,
           2500490180, 1417162592, 3978580714,  610730090,  325255158,
           2407491533, 1055922583, 3692065886, 4090616537,  326167491,
           4085233305, 2969722453, 3652458276, 2283746597, 1122234986,
            153698438, 3963478720, 1702988934, 3282802618,  486995254,
           3954255184,  643206784, 2898442177, 1561876067, 2171714320,
           3353714523, 1852240479, 3311465943, 3936017465, 3506126246,
            187682237, 2184861630, 1014617764,  632575160, 1279113960,
           3175856831, 2454753402, 1085466593,   64124646, 1379591533,
           1217101375, 2616582876, 1681054903, 2979734976, 3713477376,
           1023925369, 2241332116, 2799710965, 3289542100,  457541045,
           3394250963, 1036025802,   50390709, 3664352141, 1172119945,
           1431629985, 2978778261, 1903811186,   14623417, 1718792717,
           4034554607, 1794309220, 1156120919, 4109956536, 4213648688,
           1289890125, 1227104944, 3048900163, 1938110143, 3437073071,
           3422043195, 3967591657, 3458226637,  123675503,  779371130,
            500308599, 1670603193,  939156219, 1146229024,  757836955,
            947341990, 2325015678, 2131235396, 3980021549, 2802637798,
           1241164261, 2695616179, 2959234987,  806131530,  142482306,
           2400115798, 3226086426, 3970707668, 2061499621, 1320085716,
             57568453,  473901085,  479464733, 4011290206, 3187690206,
           1413921393, 3549162159, 2614476355, 2741496833, 2634662778,
           4181690846, 1834480108,   81313374, 3433556230,  165270435,
           2612249606, 2691848357, 2988173391, 2780259454, 3367190692,
           1237683957, 1644969050, 1948447664, 3560678650, 2791766079,
           2514319661,  266715684, 2167993587, 2420223337, 3638749523,
           1546582521,  563613841, 2217606450, 2121272794, 1342791221,
           2212917200, 2975843360, 3373518511, 4164123714, 1569112181,
            825951411, 1311236499, 3173147295,  233834319, 3423984855,
            176891089,  742869225, 1636516083, 4150973460, 3325051321,
           2781330427,  341482305,  785567544, 3505828872, 1349975937,
           2272709978, 3185457661, 1552891148, 3754754614, 2216948705,
           4200504182, 2725847895, 3319009867, 1346077197,  231273672,
            667844246, 1128800402, 2470424857, 2857709096,  939232469,
            625960350,  262896990, 2449837879, 1086596914,  526876141,
            376410310, 1945965125,  286799352, 3714930904,  605087185,
            448243262, 1200774407, 1130436969, 3242825685,   35576057,
           1654316374, 3886488657,  141637207,   64694585, 1975230700,
           2420313936, 2279935585,  703082449,  285877198, 1399629693,
           2455656193, 3473978919, 3540645678, 3556214553, 3725869341,
           3374078507, 1261526430, 3301618910, 3714779219, 1694357458,
           1003672309, 2129963260,  669089824,  952916094,  194387330,
           3438702335, 2410870937, 3822168346, 2662260186,  958314857,
            291656322,  594852624, 1571653475, 3933054030, 1444215978,
           3108222934, 2812248167, 1851787482,  583598406, 1752237696,
           1283010502, 4165771472, 1489947116, 4043644754, 2548080155,
           3528828439, 2024707308,  974264724, 1184636864, 3191666359,
              4994366, 2969723347,  203085708, 4226032208,  144035575,
           2515419987, 1748733341, 2358386124,  305081101,  134007739,
            507884251, 2196938703,  319155521, 2786035065, 2891285577,
           1182321643, 1991464298, 2518206969, 1242386622,  769489265,
           1647029574, 2077682509, 1577492020, 3569173776,  105792301,
           1448580633, 1282934123, 3488127383,  662963660, 2821302432,
           1650628199,  323225363, 1466075475, 1097647752], dtype=uint32), 175, 0, 0.0))




```python
fig, axes = plt.subplots(4, figsize=(10, 7), sharex=True)
samples = sampler.get_chain()
labels = ['A', 'B', 'C', 'nu']
for i in range(ndim):
    ax = axes[i]
    ax.plot(samples[:, :, i], "k", alpha=0.3)
    ax.set_xlim(0, len(samples))
    ax.set_ylabel(labels[i])
    ax.yaxis.set_label_coords(-0.1, 0.5)

axes[-1].set_xlabel("step number");
```


```python
flat_samples = sampler.get_chain(discard=100, thin=15, flat=True)
```


```python
fig = corner.corner(
    flat_samples, labels=labels, truths=[np.mean(A11), np.mean(B11), np.mean(C11), np.mean(Nu11)]
);
```


```python
A_emcee = flat_samples[:,0]
B_emcee = flat_samples[:,1]
nu_emcee = flat_samples[:,3]

mass_emcee = mass(A_emcee, B_emcee, nu)
mass_emcee1 = mass_emcee/5.972e+24
vari2_e = np.var(mass_emcee1)
stdd2_e = np.std(mass_emcee1)
print('Empirical Measurements:')
print('-----------------------')
print('Estimate: ' + str(np.mean(mass_emcee1)) + ' M_earth')
print('Variance: ' + str(vari2_e) + ' M_earth^2')
print('Standard Deviation: ' + str(stdd2_e) + ' M_earth')

plt.hist(mass_emcee1, bins=50)
plt.xlabel('Mass of the planet (in Earth-Mass)')
plt.ylabel('Counts')
```

a8899f6b54cbf0c7aa6ff1367437d6c7176bbb16
ipynb
Jupyter Notebook
Mass_of_exoplanet/Homework-the_mass_of_exoplanets-varying_nu.ipynb
Jayshil/Astro-data_science
8f83643197cf05981e09490352caeec3f0cde4ae
MIT
Mass_of_exoplanet/Homework-the_mass_of_exoplanets-varying_nu.ipynb
Jayshil/Astro-data_science
8f83643197cf05981e09490352caeec3f0cde4ae
MIT
Mass_of_exoplanet/Homework-the_mass_of_exoplanets-varying_nu.ipynb
Jayshil/Astro-data_science
8f83643197cf05981e09490352caeec3f0cde4ae
MIT
Qwen/Qwen-72B
1. YES
2. YES