text:
```python
%load_ext autoreload
%autoreload 2

import copy

import numpy as np
import pandas as pd
import scipy.stats as st

import re #regex

import cmdstanpy
import arviz as az

import bebi103
import bokeh_catplot

import bokeh.io
import bokeh.plotting
import bokeh.layouts
# import bokeh.models.mappers
import bokeh.palettes

import holoviews as hv
# import holoviews.operation.datashader
hv.extension('bokeh')
bebi103.hv.set_defaults()

# import panel as pn
# pn.extension()

bokeh.io.output_notebook()

%load_ext blackcellmagic
```

### Munging
Load data from Brewster, pre-tidied by Manuel, and drop the spurious column that was the index in csv.
See `code/exploratory/fish_munging.ipynb` for details. TL;DR: don't use the regulated csv, the one below has all the FISH data. mRNA_cell is the data we want, not spots_totals.

```python
df_fish = pd.read_csv("../../data/jones_brewster_2014.csv")
del df_fish['Unnamed: 0']
df_fish.head()
```

#### Separating out regulated data
The regulated datasets' labels start with O1, O2, or O3. Everything else doesn't. Use regex to parse.

```python
raw_expt_labels = df_fish['experiment'].unique()
raw_expt_labels.sort()

# put all strings that start w/ 'O' in one list
regulated_labels = [label for label in raw_expt_labels if re.match('^O', label)]
```

This leaves behind some data with insufficient metadata, e.g., what good does the aTc concentration do me if I don't know what promoter it was for?

Now that we've got labels we want, let's slice dataframes accordingly. O1, 2ng/mL will be our test sample for SBC.

```python
df_reg = df_fish[df_fish['experiment'].isin(regulated_labels)]
df_UV5 = df_fish[df_fish["experiment"] == "UV5"]
df_rep_test = df_reg[df_reg['experiment'] == 'O1_2ngmL']
```

## Analyzing simple repression with Stan

Stan model borrows from JB's tutorial 7a, 2018, and from JB's finch beak tutorial for bebi103b 2020 TAs.

The model here is the same negative binomial model we used for the constitutive case, except with the burst rate multiplied by the fold-change in mean. Can we successfully infer the Bohr parameter?

As in the constitutive case, Stan parametrizes the negative binomial with $\alpha$ and $\zeta$, where $\alpha$ is the burst frequency (dimensionless, nondimensionalized by mRNA lifetime) and $\zeta = 1/b$ where $b$ is the mean burst size.

What is our prior range for fold-change? Consider
\begin{align}
fc = \left( 1+\frac{R}{N_{NS}} e^{-\beta\Delta\epsilon} \right)^{-1}
   = \left( 1+
        \exp\left(\ln\left(\frac{R}{N_{NS}}\right)-\beta\Delta\epsilon\right)
     \right)^{-1}
   = \left( 1+ e^{-\beta\Delta F} \right)^{-1}
\end{align}
$R/N_{NS}$ should definitely remain between $10^{-7}$ and $10^{-2}$, and we expect/know the $\Delta\epsilon$ to range between about $-18$ and $-9$. So we can say $\beta\Delta F$ is unlikely to escape a range between $-14$ and $+7$, corresponding to fold-changes of $10^{-6}$ and $(1-10^{-3})$, which certainly already exceed our detection limits anyways. A prior such as
\begin{align}
\beta\Delta F \sim \text{Normal}(\mu=-3,\sigma=5)
\end{align}
corresponds to this story.

Then our target likelihood is
\begin{align}
p_m \sim \text{NBinom}(\alpha, \zeta)
\end{align}
for the UV5 data, and simultaneously
\begin{align}
p_m \sim \text{NBinom}\left(\frac{\alpha}{1+ e^{-\beta\Delta F}}, \zeta\right)
\end{align}
for the regulated data.

### Prior predictive checks

Note that my priors on $\alpha$ and $b$ are tighter than what I used for the constitutive analysis. That is by design, since we now know approximately what they should be for UV5. They are still weaker than the posterior was for constitutive UV5, but this should help prevent the repressed data from pulling them to crazy values unless the data really warrants it.

```python
sm_prior_predictive = cmdstanpy.CmdStanModel(
    stan_file="stan/simple_rep_means_prior_predictive_v01.stan"
)
# print(sm_prior_predictive.code())
```

```python
data_prior_pred = dict(
    N_cells_uv5=len(df_UV5),
    N_cells_rep=len(df_rep_test),
    log_alpha_loc=1.0,
    log_alpha_scale=1.0,
    log_b_loc=1.0,
    log_b_scale=0.5,
    bohr_loc=-3.0,
    bohr_scale=3.0
)
```

```python
prior_pred_samples = sm_prior_predictive.sample(
    data=data_prior_pred,
    fixed_param=True,
    sampling_iters=1000,
    output_dir="./stan/stan_samples",
)
```

```python
# Convert to ArviZ InferenceData
prior_pred_samples = az.from_cmdstanpy(
    posterior=prior_pred_samples,
    prior=prior_pred_samples,
    prior_predictive=['mRNA_counts_uv5', 'mRNA_counts_rep']
)
```

```python
p1 = bebi103.viz.predictive_ecdf(
    prior_pred_samples.prior_predictive['mRNA_counts_uv5'],
    frame_height=250,
    frame_width=350,
    discrete=True,
    percentiles=(95, 90, 75, 50),
    x_axis_label='mRNA counts (UV5)',
    x_axis_type='log'
)
p2 = bebi103.viz.predictive_ecdf(
    prior_pred_samples.prior_predictive['mRNA_counts_rep'],
    frame_height=250,
    frame_width=350,
    discrete=True,
    percentiles=(95, 90, 75, 50),
    x_axis_label='mRNA counts (repressed)',
    x_axis_type='log'
)
bokeh.io.show(p1)
bokeh.io.show(p2)
```

### Simulation-based calibration

Next up: simulation-based calibration (SBC). Quoting JB, this checks "that the sampler can effectively sample the entire space of parameters covered by the prior." We'll go ahead and set up the data for the posterior, even though we won't be sampling the posterior right now. Then we can set up the model.

```python
data_rep_test = copy.deepcopy(data_prior_pred)
data_rep_test["N_cells_rep"] = len(df_rep_test)
data_rep_test["mRNA_counts_uv5"] = df_UV5["mRNA_cell"].values.astype(int)
data_rep_test["mRNA_counts_rep"] = df_rep_test["mRNA_cell"].values.astype(int)
data_rep_test["ppc"] = 0

sm = cmdstanpy.CmdStanModel(stan_file="stan/simple_rep_means_v01.stan")
# print(sm.code())
```

Now just run it. If for some reason you're running on a laptop (why?) rather than Delbruck, don't forget to reduce # of cores.

```python
try:
    sbc_output = pd.read_csv("stan/sbc_simple_rep_means_v01.csv")
except:
    sbc_output = bebi103.stan.sbc(
        prior_predictive_model=sm_prior_predictive,
        posterior_model=sm,
        prior_predictive_model_data=data_prior_pred,
        posterior_model_data=data_rep_test,
        measured_data=["mRNA_counts_uv5", "mRNA_counts_rep"],
        parameters=["alpha", "b", "bohr"],
        sampling_kwargs={'thin': 10},
        cores=45,
        N=400,
        progress_bar=True,
    )

    sbc_output.to_csv("stan/sbc_simple_rep_means_v01.csv", index=False)
```

Plot ECDFs of the rank statistics, which should be ~ uniform. Color according to warning code, which will also let us know if there were any issues with divergences, R-hat, EBFMI, effective number of steps, or tree depth.

```python
plots = [
    bokeh_catplot.ecdf(
        data=sbc_output.loc[sbc_output["parameter"] == param, :],
        val="rank_statistic",
        cats="warning_code",
        kind="colored",
        frame_width=400,
        frame_height=150,
        title=param,
        conf_int=True,
    )
    for param in sbc_output["parameter"].unique()
]

bokeh.io.show(bokeh.layouts.gridplot(plots, ncols=1))
```

Warning code 2 means R-hat failure. Let's look at the R-hat values.

```python
bokeh.io.show(
    bokeh_catplot.ecdf(data=sbc_output, val="Rhat", cats="parameter")
)
```

If there are a substantial fraction of Rhats above the brightline 1.01, JB suggests that more samples from the posterior might be wise.

JB also has a nice helper function to plot the difference from uniform of rank ECDFs, with a confidence interval for the uniform distribution to guide the eye. Let's take a look.

```python
bokeh.io.show(bebi103.viz.sbc_rank_ecdf(sbc_output))
```

Finally let's check shrinkage and z-scores.

```python
hv.Points(
    data=sbc_output,
    kdims=['shrinkage', 'z_score'],
    vdims=['parameter', 'ground_truth', 'mean', 'sd'],
).opts(
    color='parameter',
    alpha=0.3,
    xlim=(0, 1.05),
    tools=['hover']
)
```

This actually looks better than I expected. Most of the cases of poor shrinkage on Bohr parameter occur when expected: very strong or very weak repression (e.g., Bohr $\lesssim -6~k_BT$ or  $\gtrsim 3~k_BT$), so there is not much signal to inform the inference of Bohr either way. I expected there to be more |z-scores| $\gtrsim 2$ on $\alpha$ and $b$, to reflect the pull that we see in the posteriors below, but apparently that's not the case after all.

### Sampling the Posterior
We already finished building the model in order to do SBC. Now we just run it.

```python
# We do want posterior predictive checks this time
data_rep_test["ppc"] = 1

posterior_samples = sm.sample(data=data_rep_test, cores=5)
posterior_samples = az.from_cmdstanpy(
    posterior_samples, posterior_predictive=["mRNA_counts_uv5_ppc", "mRNA_counts_rep_ppc"]
)
```

```python
bebi103.stan.check_all_diagnostics(posterior_samples)
```

Let's visualize the posterior.

```python
bokeh.io.show(
    bebi103.viz.corner(
        posterior_samples,
        pars=["alpha", "b", "bohr"],
        alpha=0.1,
        xtick_label_orientation=np.pi / 4,
    )
)
```

That looks quite reasonable. The transcripts per burst & burst frequency are both comparable to what we would have inferred from Manuel's MCMC, and also what we inferred from UV5 data alone.

### Posterior predictive checks

Even though the posterior looks ok, the model could still be wrong: it is identifiable, but is it consistent with the data? Posterior predictive checks address this by asking whether the model could plausibly generate the observed data. (Function borrowed from JB's finch beak tutorial for bebi103b 2020 TAs.)

```python
def ppc_ecdf_pair(posterior_samples, ppc_var, df, percentiles=(80, 60, 40, 20),
                 x_axis_label="mRNA counts per cells", frame_width=200, frame_height=200):
    """Plot posterior predictive ECDFs."""
    n_samples = (
        posterior_samples.posterior_predictive.dims["chain"]
        * posterior_samples.posterior_predictive.dims["draw"]
    )

    p1 = bebi103.viz.predictive_ecdf(
        posterior_samples.posterior_predictive[ppc_var].values.reshape(
            (n_samples, len(df))
        ),
        data=df["mRNA_cell"],
        percentiles=percentiles,
        discrete=True,
        x_axis_label=x_axis_label,
        frame_width=frame_width,
        frame_height=frame_height
    )

    p2 = bebi103.viz.predictive_ecdf(
        posterior_samples.posterior_predictive[ppc_var].values.reshape(
            (n_samples, len(df))
        ),
        data=df["mRNA_cell"],
        percentiles=percentiles,
        discrete=True,
        x_axis_label=x_axis_label,
        frame_width=frame_width,
        frame_height=frame_height,
        diff=True,
    )
    p1.x_range = p2.x_range

    return [p1, p2]
```

```python
ppc_plots = (ppc_ecdf_pair(posterior_samples, "mRNA_counts_uv5_ppc", df_UV5),
             ppc_ecdf_pair(posterior_samples, "mRNA_counts_rep_ppc", df_rep_test))

# flatten list of lists w/ list comp, then plot
ppc_plots = [subplot for sublist in ppc_plots for subplot in sublist]
bokeh.io.show(bokeh.layouts.gridplot(ppc_plots, ncols=2))
```

Depending on which lac operator and aTc concentration we analyze, I'd draw rather different conclusions. Let's just analyze everything.

### Sampling all the data
Let's repeat for all the operators and repressor copy numbers! Since we have so many, do separate loops to generate the samples and generate viz (so we can tweak viz without pointlessly rerunning the sampling). Put all the samples in a dict for easy access; I'm not sure if arviZ will keep them in ram, but model is small enough that even if so, it'll be fine.

```python
all_samples = {}
for trial in df_reg["experiment"].unique():
    temp_df = df_reg[df_reg["experiment"] == trial]
    data = copy.deepcopy(data_prior_pred)
    data["N_cells_uv5"] = len(df_UV5)
    data["N_cells_rep"] = len(temp_df)
    data["mRNA_counts_uv5"] = df_UV5["mRNA_cell"].values.astype(int)
    data["mRNA_counts_rep"] = temp_df["mRNA_cell"].values.astype(int)
    data["ppc"] = 1

    # disable_logging suppresses cmdstanpy output, uncomment if you want to watch it sample
    with bebi103.stan.disable_logging():
        posterior_samples = sm.sample(data=data, cores=6)

    all_samples[trial] = az.from_cmdstanpy(
        posterior_samples,
        posterior_predictive=["mRNA_counts_uv5_ppc", "mRNA_counts_rep_ppc"],
    )
```

It's always wise to check diagnostics.

```python
# for trial in all_samples:
#     print(trial)
#     bebi103.stan.check_all_diagnostics(all_samples[trial])
```

There's some run-to-run variation. Sometimes one or two promoters give Rhat warnings a bit over 1.01, but generally this looks ok.

Now that we've sampled we can plot. (I can't figure out how to add titles to the bokeh layouts that `viz.corner` returns, so I'm doing the poor man's version for now.)

```python
for trial in all_samples:
    print(trial)
    # plot posterior as corner plot
    bokeh.io.show(
        bebi103.viz.corner(
            all_samples[trial],
            pars=["alpha", "b", "bohr"],
            alpha=0.1,
            xtick_label_orientation=np.pi / 4,
        )
    )
    # setup post pred ecdfs
    ppc_plots = (
        ppc_ecdf_pair(all_samples[trial], "mRNA_counts_uv5_ppc", df_UV5),
        ppc_ecdf_pair(
            all_samples[trial],
            "mRNA_counts_rep_ppc",
            df_reg[df_reg["experiment"] == trial],
        ),
    )

    # flatten list of lists to prepare for gridplot
    ppc_plots = [subplot for sublist in ppc_plots for subplot in sublist]
    bokeh.io.show(bokeh.layouts.gridplot(ppc_plots, ncols=2))
```

First off, obviously this model is not "true." But it's not a bad zeroth-order model for simple repression: it provides an interesting way to link the thermodynamic model (Bohr parameter) to observables beyond means.

And once again this model performs better than I'd expect it to, and it's failure mode is very enlightening. It is _extremely_ interesting that for the very weakly repressed and the very highly repressed conditions, the joint posteriors on $\alpha$ and $b$ (i.e., after marginalizing away Bohr) are _very_ similar to the posterior we got from the UV5 data alone. In these limits, approximating the distribution as still being negative binomial with a rescaled burst rate seems alright. Another statement of this: in these limits, there is still only one timescale in sight: the burst rate, appropriately scaled.

In the intermediate regime, this is not true: the repressed mRNA distributions are substantially more disperse than constitutive UV5. The model can't fit both dists simultaneously, so it contorts itself to try and split the difference. `Oid0p5ngmL` is probably the most extreme example. $\alpha$ gets pulled down and $b$ gets pulled up to increase the variance, resulting in posterior predictive distributions that are too wide for UV5 and still too narrow for the repressed data.

This is actually _a good sign_: the addition of repressor leaves a nontrivial imprint on the mRNA distribution. So a kinetic model with repressor dynamics may actually be identifiable. The question is whether, even if we can infer 2 rates and label them $k_R^+$ and $k_R^-$, do they actually correspond to the microscopic repressor kinetics, or are they just related to cell cycle or mRNA lifetimes or TetR partitioning imprecision or or...? But we'll cross that bridge when we get there...

#### Inferring repressor copy number
If we take the binding energies of each operator as knowns, we can extract $R/N_{NS}$ for each op/aTc pair and see if they make sense. First collect the approximate MAP of the Bohr parameter (in $kT$) for each trial, tabulated below.

Inferred Bohr parameters:

|     | 0.5 ng/ml | 1 ng/ml | 2 ng/ml | 10 ng/ml |
| --- | --------- | ------- | ------- | -------- |
| Oid | 1.2       | -2.75   | -4.8    | -4.4     |
| O1  | 2.5       | -1.8    | -3.45   | -3.6     |
| O2  | 1.8       | -0.5    | -1.8    | -2.8     |
| O3  | ?         | > 3     | > 2.5   | 2 ?      |

Then using
\begin{align}
\Delta F = \Delta \epsilon - \ln(R/N_{NS})
\end{align}
and
\begin{align}
\Delta \epsilon_{Oid} = -17.7 \\
\Delta \epsilon_{O1}  = -15.3 \\
\Delta \epsilon_{O2}  \in (-13.9, -13.6) \\
\Delta \epsilon_{O3}  \in (-9.7, -9.4) \\
\end{align}
we get estimates for $R$.

Inferred $R$:

|     | 0.5 ng/ml | 1 ng/ml | 2 ng/ml | 10 ng/ml |
| --- | --------- | ------- | ------- | -------- |
| Oid | 0.03      | 1.5     | 11      | -        |
| O1  | 0.1       | 6       | 30-35   | -        |
| O2  | 0.5-1     | 6-8     | 25-35   | 70-100   |
| O3  | -         |  -      | > 25    | > 30?    |

We have left blank values where the Bohr parameter was clearly inferred poorly (either very strong or very weak repression), which would result in obviously inaccurate estimates of $R$. Even the values listed should be taken with a grain of salt, but this gives us an order of magnitude sense of things.

Strangely, though, Brewster and Daniel were using HG203 as a base just as I am, with _tetR_ integrated at the _gspI_ locus and _lacI_ at _ybcN_. The only difference is the reporter at _galK_ and the mNeonGreen fusion. So that means our aTc induction levels ought to be approximately apples to apples, no? Why do I seem to get $\sim1$ order of magnitude stronger induction?? Or maybe my estimates here are just that sloppy?

#### Next step

So we can fit each new repressor copy number & operator combo with it's own _de novo_ Bohr parameter & that is an easily identifiable model. Good.

Can we be more principled? For the next iteration of the model, let's define a binding energy for each operator and a single fit parameter that globally converts from aTc to lacI copy number (assume they're linearly proportional). Fit all the data at once. Now there's only one global $\alpha$ and $b$; how will they adapt? And can we recover, at least approximately, the known operator binding energies.

For coding the full true likelihood, JB suggests referring to a JavaScript library of special functions, by paulmasson.


text:
# Almgren and Chriss Model For Optimal Execution of Portfolio Transactions

### Introduction

We consider the execution of portfolio transactions with the aim of minimizing a combination of risk and transaction costs arising from permanent and temporary market impact. As an example, assume that you have a certain number of stocks that you want to sell within a given time frame. If you place this sell order directly to the market as it is, transaction costs may rise due to temporary market impact. On the other hand, if you split up into pieces in time, cost may rise due to volatility in the stock price.

[Almgren and Chriss](https://cims.nyu.edu/~almgren/papers/optliq.pdf) provided a solution to this problem by assuming the permanent and temporary market impact functions are linear functions of the rate of trading, and that stock prices follow a discrete arithmetic random walk.

In this notebook, we will take a look at the model used by Almgren and Chriss to solve the optimal liquidation problem. We will start by stating the formal definitions of *trading trajectory*, *trading list*, and *trading strategy* for liquidating a single stock.

### Trading Trajectory, Trading List, and Trading Strategy

We define trading trajectory, trading list, and trading strategy just as Almgren and Chriss did in their [paper](https://cims.nyu.edu/~almgren/papers/optliq.pdf). Suppose we hold $X$ shares of a stock that we want to liquidate before time $T$. Divide $T$ into $N$ intervals of length $\tau=\frac{T}{N}$ and define:

- $t_k = k\tau$ to be discrete times, where $k = 0,..,N$.

- A **trading trajectory** to be the list $(x_0,..,x_N)$, where $x_k$ is the number of shares we plan to hold at time $t_k$. We require that our initial position $x_0 = X$, and that at liquidation time $T$, $x_N = 0$.

- A **trading list** to be $(n_1,..,n_N)$, $n_k =  x_{k-1} - x_k$ as the number of shares that we will sell between times $t_{k-1}$ and $t_k$.

- A **trading strategy** as a rule for determining $n_k$ from the information available at time $t_{k-1}$.

Below, we can see a visual example of a trading trajectory, for $N = 12$.

## Price Dynamics

We will assume that the stock price evolves according to a discrete arithmetic random walk:

\begin{equation}
S_k = S_{k-1} + \sigma \tau^{1/2} \xi_k
\end{equation}

for $k = 1,..,N$ and where:

\begin{equation}
S_k = \text{ stock price at time $k$}\hspace{21.6cm}\\
\sigma = \text{ standard deviation of the fluctuations in stock price}\hspace{16.3cm}\\
\tau = \text{ length of discrete time interval}\hspace{20.2cm}\\
\xi_k = \text{ draws from independent random variables}\hspace{17.8cm}
\end{equation}

We will denote the initial stock price as $S_0$. The role of $\xi_k$ is to simulate random price fluctuations using random numbers drawn from a Normal Gaussian distribution with zero mean and unit variance. The code below shows us what this price model looks like, for an initial stock price of $S_0 =$ \$50$ dollars, a standard deviation of price fluctuations of $\sigma = 0.379$, and a discrete time interval of $\tau = 1$.

```python
%matplotlib inline

import matplotlib.pyplot as plt

import utils

# Set the default figure size
plt.rcParams['figure.figsize'] = [17.0, 7.0]

# Set the number of days to follow the stock price
n_days = 30

# Plot the stock price as a function of time
utils.plot_price_model(seed = 3, num_days = n_days)
```

## Market Impact

As we learned previously the price of a stock is affected by market impact that occurs every time we sell a stock. In their model, Almgren and Chriss distinguish between two types of market impact, permanent and temporary market impact. We will now add these two factors into our price model.

### Permanent Impact

Permanent market impact refers to changes in the equilibrium price of a stock as a direct function of our trading. Permanent market impact is called *permanent* because its effect persists for the entire liquidation period, $T$. We will denote the permanent price impact as $g(v)$, and will add it to our price model:

\begin{equation}
S_k = S_{k-1} + \sigma \tau^{1/2} \xi_k - \tau g\left(\frac{n_k}{\tau}\right)
\end{equation}

Here, we assumed the permanent impact function, $g(v)$, is a linear function of the trading rate, $v = n_k / \tau$. We will take $g(v)$ to have the form:

\begin{equation}
g(v) = \gamma \left(\frac{n_k}{\tau}\right)
\end{equation}

where $\gamma$ is a constant and has units of \\$ per ${share}^{2}$. Replacing this in the above equation we get:

\begin{equation}
S_k = S_{k-1} + \sigma \tau^{1/2} \xi_k - \gamma n_k
\end{equation}

With this form, we can see that for each $n$ shares that we sell, we will depress the stock price permanently by $n\gamma$, regardless of the time we take to sell the stocks.

### Temporary Impact

Temporary market impact refers to temporary imbalances in supply and demand caused by our trading. This leads to temporary price movements away from equilibrium. Temporary market impact is called *temporary* because its effect
dissipates by the next trading period. We will denote the temporary price impact as $h(v)$. Given this, the actual stock price at time $k$ is given by:

\begin{equation}
\tilde{S_k} = S_{k-1} - h\left(\frac{n_k}{\tau}\right)
\end{equation}

Where, we have again assumed the temporary impact function, $h(v)$, is a linear function of the trading rate, $v = n_k / \tau$. We will take $h(v)$ to have the form:

\begin{equation}
h(v) = \epsilon \mbox{ sign}(n_k) + \eta \left(\frac{n_k}{\tau}\right)
\end{equation}

where $\epsilon$ and $\eta$ are constants with units (\$/share) and (\$ time/share${}^2$), respectively. It is important to note that $h(v)$ does not affect the price $S_k$.

## Capture

We define the **Capture** to be the total profits resulting from trading along a particular trading trajectory, upon completion of all trades. We can compute the capture via:

\begin{equation}
\sum\limits_{k=1}^{N} n_k \tilde{S_k} = X S_0 + \sum\limits_{k=1}^{N} \left(\sigma \tau^{1/2} \xi_k - \tau g\left(\frac{n_k}{\tau}\right)\right) x_k - \sum\limits_{k=1}^{N} n_k h\left(\frac{n_k}{\tau}\right)
\end{equation}

As we can see this is the sum of the product of the number of shares $n_k$ that we sell in each time interval, times the effective price per share $\tilde{S_k}$ received on that sale.

## Implementation Shortfall

We define the **Implementation Shortfall** as the total cost of trading and is given by:

\begin{equation}
I_s = X S_0 - \sum_{k = 1}^N n_k \tilde{S_k}
\end{equation}

This is what we seek to minimize when determining the best trading strategy!

Note that since $\xi_k$ is random, so is the implementation shortfall. Therefore, we have to frame the minimization problem in terms of the expectation value of the shortfall and its corresponding variance. We'll refer to $E(x)$ as the expected shortfall and $V(x)$ as the variance of the shortfall. Simplifying the above equation for $I_s$, is easy to see that:

\begin{equation}
E(x) = \sum\limits_{k=1}^{N} \tau x_k g\left(\frac{n_k}{\tau}\right) + \sum\limits_{k=1}^{N} n_k h\left(\frac{n_k}{\tau}\right)
\end{equation}

and

\begin{equation}
V(x) = \sigma^2 \sum\limits_{k=1}^{N} \tau {x_k}^2
\end{equation}

The units of $E(x)$ are dollars and the units of $V(x)$ are dollars squared. So now, we can reframe our minimization problem in terms of $E(x)$ and $V(x)$.
For a given level of variance of shortfall, $V(x)$, we seek to minimize the expectation of shortfall, $E(x)$. In the next section we will see how to solve this problem.

## Utility Function

Our goal now is to find the strategy that has the minimum expected shortfall $E(x)$ for a given maximum level of variance $V(x) \ge 0$. This constrained optimization problem can be solved by introducing a Lagrange multiplier $\lambda$. Therefore, our problem reduces to finding the trading strategy that minimizes the **Utility Function** $U(x)$:

\begin{equation}
U(x) = E(x) + \lambda V(x)
\end{equation}

The parameter $\lambda$ is referred to as **trader’s risk aversion** and controls how much we penalize the variance relative to the expected shortfall.

The intuition of this utility function can be thought of as follows. Consider a stock which exhibits high price volatility and thus a high risk of price movement away from the equilibrium price. A risk averse trader would prefer to trade a large portion of the volume immediately, causing a known price impact, rather than risk trading in small increments at successively adverse prices. Alternatively, if the price is expected to be stable over the liquidation period, the trader would rather split the trade into smaller sizes to avoid price impact. This trade-off between speed of execution and risk of price movement is ultimately what governs the structure of the resulting trade list.

# Optimal Trading Strategy

Almgren and Chriss solved the above problem and showed that for each value
of risk aversion there is a uniquely determined optimal execution strategy. The details of their derivation is discussed in their [paper](https://cims.nyu.edu/~almgren/papers/optliq.pdf). Here, we will just state the general solution.

The optimal trajectory is given by:

\begin{equation}
x_j = \frac{\sinh \left( \kappa \left( T-t_j\right)\right)}{ \sinh (\kappa T)}X, \hspace{1cm}\text{ for } j=0,...,N
\end{equation}

and the associated trading list:

\begin{equation}
n_j = \frac{2 \sinh \left(\frac{1}{2} \kappa \tau \right)}{ \sinh \left(\kappa T\right) } \cosh \left(\kappa \left(T - t_{j-\frac{1}{2}}\right)\right) X, \hspace{1cm}\text{ for } j=1,...,N
\end{equation}

where $t_{j-1/2} = (j-\frac{1}{2}) \tau$. The expected shortfall and variance of the optimal trading strategy are given by:

In the above equations $\kappa$ is given by:

\begin{align*}
&\kappa = \frac{1}{\tau}\cosh^{-1}\left(\frac{\tau^2}{2}\tilde{\kappa}^2 + 1\right)
\end{align*}

where:

\begin{align*}
&\tilde{\kappa}^2 = \frac{\lambda \sigma^2}{\tilde{\eta}} = \frac{\lambda \sigma^2}{\eta \left(1-\frac{\gamma \tau}{2 \eta}\right)}
\end{align*}

```python

```

ext: ipynb

max_stars_repo_path: finance/AlmgrenChriss_notebook/Almgren and Chriss Model.ipynb


text:
```python
import random
from sympy import *
```

```python
Matrix??
```

```python
B = Matrix(3,3, lambda i,j: S(random.randrange(-1,4))/random.randrange(1,4) )
B
```

    Matrix([
    [   0,   1, 1/3],
    [-1/3, 1/3, 2/3],
    [ 2/3,   1,   2]])

```python
B.det()
```

    25/27

```python

```

```python

```

```python
beta = symbols('beta')

```

```python
A
```

```python
A.det()
```

```python
# if alpha = 2/3 then det(A)=0  => A is not invertible
# otherwise A is invertible.

```

```python
S(random.randrange(1,4))/random.randrange(2,4)
```

```python
A = Matrix([[3,1],[2,S(3)/2]])
```

```python
A

```

```python

```

max_stars_repo_name: ChidinmaKO/noBSLAnotebooks

max_issues_repo_name: ChidinmaKO/noBSLAnotebooks

max_forks_repo_name: ChidinmaKO/noBSLAnotebooks


text:
# Lecture 28

## Sum of a random number of random variables; Inequalities (Cauchy-Schwarz, Jensen, Markov, Chebyshev)
----

```python
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import beta
%matplotlib inline

plt.xkcd()
```

    <contextlib._GeneratorContextManager at 0x15cd2c19550>

## Sum of Random Number of Random Variables

Let's say there's a store. At this store, there is a random number of customers. Each customer chooses to spend a random amount of money at the store. What can we say about the expected total amount spent at the store? What about its variance?

* $N$ is the number of customers (maybe Poisson? it doesn't matter at this point)
* $X_j$ be the amount customer $j$ spends
* $X_j$ has mean $\mu$ and variance $\sigma^2$
* assume that $N, \, X_1, \, N_2, \dots$ are independent

Find the mean and variance of $X = \sum_{j=1}^{N} X_j$.

Be careful, because $N$ is itself a random variable, and so we are summing up a _random number of random variables_. At this point, $N$ could be any distribution, perhaps Poisson, but we can put that question aside for now. Let's see if we can find some statistics about $X$ in terms of random variable $N$.

### $\mathbb{E}(X)$

Now, if we blindly tried to apply linearity, we might have:

\begin{align}
  \mathbb{E}(X) &=? \, N \, \mu &\quad \text{... no, this cannot be right}
\end{align}

But this is a category error, since the right-hand side of that equation is a _value_, while the left-hand side is a _random variable_. However, it does give us a clue: all we need to do is find a way to treat $N$ with a _value_, which we can do by thinking conditionally.

There are two ways we can use to find $\mathbb{E}(X)$:

\begin{align}
  &\text{...conditioning on } N \text{, the long way} \\
  \mathbb{E}(X) &= \sum_{n=0}^{\infty} \mathbb{E}(X|N=n) \, P(N=n) &\quad \text{conditioning on }N \\
  &= \sum_{n=0}^{\infty} \mu \, n \, P(N=n) &\quad \text{by linearity, independence} \\
  &= \mu \, \sum_{n=0}^{\infty} n \, P(N=n) &\quad \text{but by definition of }\mathbb{E} \\
  &= \boxed{\mu \, \mathbb{E}(N)} \\\\
  \\
  & \text{... or using Adam's Law} \\
  \mathbb{E}(X) &= \mathbb{E}\left( \mathbb{E}(X|N) \right) &\quad \text{Adam's Law} \\
  &= \mathbb{E}(\mu \, N) &\quad \text{treating } N \text{ as constant, and by linearity} \\
  &= \boxed{\mu \, \mathbb{E}(N)}
\end{align}

Note that unlike the 2-envelope paradox covered in lectures 25 and 26, we _can_ forget about the condition $N=n$ _because of the independence of $X_j$_. That fact allows us to use linearity in the second step of finding $\mathbb{E}(X)$ the long way by conditioning on $N$.

### $Var(X)$

Here, we will just apply EVvE's Law:

\begin{align}
  Var(X) &= \mathbb{E}\left( Var(X|N) \right) + Var\left( \mathbb{E}(X|N) \right) \quad \text{and working on each part} \\\\
  \mathbb{E}\left( Var(X|N) \right) &= \mathbb{E}(N \sigma^2) \quad \text{treating } N \text{ as fixed, independence} \\
  &= \sigma^2 \, \mathbb{E}(N) \\
  \\
  Var\left( \mathbb{E}(X|N) \right) &= Var(\mu \, N) \quad \text{} \\
  &= \mu^2 \, Var(N) \\\\
  \Rightarrow Var(X) &= \boxed{ \sigma^2 \, \mathbb{E}(N) + \mu^2 \, Var(N) }
\end{align}

Idiot-checking that result, we confirm that

* $\mathbb{E}(X) = \mu \, \mathbb{E}(N)$ means that the average total sales is the average per-customer spending amount times the average number of customers
* $\sigma^2$ is in terms of money-squared; $\mu$ is in terms of money, so $Var(X)$ is in terms of money-squared
* both $\mathbb{E}(N)$ and $Var(N)$ are dimensionless values ("people" are not units in this case)
* if we knew the MGF of each $X_j$, then we could condition on $N$ to treat it as fixed in order to get the MGF of $X$ (multiply by $n$)

## Statistical Inequality

In the field of law, it is a lot easier to stand in court to use statistical inequalities rather than approximations.

Basically, approximations could be painted to be subjective. If you use a Poisson approximation and had to explain yourself in court, you might be hard-pressed to justify why you think your approximation is good. "Good" is subject to interpretation; we would claim that an approximation is good if it is close to the truth. But a lawyer might then grill on that, asking you "Is there some accepted standard to closeness? Do you know how close your approximation is?"

Now if you knew that, then you would already *know* the actual distribution and you wouldn't be using an approximation in the first place!

...

That sort of argument probably would not fly very well in front of the judge.

But inequalities are pretty much unassailable, as they provide a _definite fact about something that is random_.

### (1) Cauchy-Schwarz Inequality: $\lvert \mathbb{E}(XY) \rvert \le \sqrt{\mathbb{E}(X^2)(Y^2)}$

Recall that $\mathbb{E}(XY)$ is an analog of the dot-product for random variables $X,Y$.

Recall the definition of a dot-product in linear algebra.

\begin{align}
  \textbf{a} \cdot \textbf{b} &= \lVert \mathbf{a} \rVert \lVert \mathbf{b} \rVert \mathbf{cos}\theta \\
  \\
  \mathbb{E}(XY) &\le \sqrt{\mathbb{E}(X^2)(Y^2)}
\end{align}

Note that $X,Y$ are uncorrelated, then by definition $\mathbb{E}(XY) = \mathbb{E}(X)\mathbb{E}(Y)$, and so using the Cauchy-Schwarz Inequality is really unnecessary.

The Cauchy-Schwarz Inequality really comes in handy when the random variables $X,Y$ are correlated. In such a case where $X,Y$ are correlated and you are interested in $\mathbb{E}(XY)$, then without using a 2D LOTUS or calculating the distribution of $XY$, we can still find an upper bound to $\mathbb{E}(XY)$.

#### Statistical Interpretation

The interpretation is easiest to see if we let $X,Y$ have mean 0.

Then:
\begin{align}
  \lvert Corr(X,Y) \rvert &= \lvert \frac{Cov(X,Y)}{\sqrt{Var(X)} \sqrt{Var(Y)}} \rvert \\
  &= \lvert \frac{\mathbb{E}(XY) - \mathbb{E}(X)\mathbb{E}(Y)}{\sqrt{\mathbb{E}(X^2)-(\mathbb{E}(X))^2}\sqrt{\mathbb{E}(Y^2)-(\mathbb{E}(Y))^2}} \rvert \\
  &= \lvert \frac{\mathbb{E}(XY)}{\left(\mathbb{E}(X^2) \mathbb{E}(Y^2)\right)^{\frac{1}{2}}} \rvert \le 1
\end{align}

... and the left-hand side of the above equation is *exactly* the Cauchy-Schwarz Inequality. Note that this would hold for any value of mean for $X,Y$.

So for statistics, the Cauchy-Schwarz Inequality really just means that the correlation of two random variables will be between -1 and 1. For providing upper bounds to $\mathbb{E}(XY)$, the Cauchy-Schwarz Inequality's strengths are simplicity and generality.

### (2) Jensen's Inequality: if $g$ is a convex function, then $\mathbb{E} g(X) \ge g(\mathbb{E}X)$

\begin{align}
  \mathbb{E}g(x) &\le g(\mathbb{E}X) &\quad \text{if }g(x) \text{ is concave} \\\\
  \mathbb{E}h(x) &\ge h(\mathbb{E}X) &\quad \text{if }h(x) \text{ is convex}
\end{align}

This tells you when $g$ is a convex function which way the inequality goes.

Now, a convex function is such that $g^{\prime\prime}(x) \ge 0$ (if $g^{\prime\prime}(x)$ exists). Another way of describing it is that if you pick any two points on $g(x)$, the line segment they make will lay above $g(x)$.

```python
x = np.linspace(-10, 10, 1000)
y = x**2
x2 = [-4,7]
y2 = [16,49]

_, ax = plt.subplots(figsize=(8,6))
ax.axhline(y=0, color='k')
ax.axvline(x=0, color='k')

ax.set_xticks([])
ax.set_yticks([])

plt.text(2.0, 80.0, '$g(x)=x^2$,\n$g^{\prime\prime}(x)=2$', color='xkcd:azure')
plt.text(-5.5, 10.0, 'A')
plt.text(7.5, 42.5, 'B')

ax.plot(x, y, 'xkcd:azure', lw=4, alpha=0.6)
ax.plot(x2,y2,'xkcd:orange', linestyle='--', lw=4, alpha=0.6)

ax.title.set_text('Convexity: $\overline{AB}$ lies above $g(x)$ for all points $A,B$')
```

```python
x = np.linspace(0.001, 5, 1000)
y = 1/x
x2 = np.linspace(0.001,5,1000)
y2 = np.log(x2)

_, ax = plt.subplots(figsize=(8,6))
ax.axhline(y=0, color='k')
ax.axvline(x=0, color='k')

ax.set_xticks([])
ax.set_yticks([])
ax.set_xlim([-2,5])
ax.set_ylim([-2,5])
ax.text(0.85, 2.2, (r'$\mathbb{E}(\frac{1}{X}) \geq \mathbb{E}(X)$' '\n$X$ is positive, $g(X)$ is convex'), color='xkcd:azure')
ax.text(0.85, -1.2, '$\mathbb{E}(\ln X) \leq \ln \mathbb{E}(X)$\n$h(X)$ is concave, inequality flips', color='xkcd:orange')
ax.plot(x, y, 'xkcd:azure', lw=4, alpha=0.6)
ax.plot(x2, y2, 'xkcd:orange', lw=4, alpha=0.6)

ax.title.set_text("Examples: Jensen's Inequality")
```

#### Proof

We will try a graphical approach to prove Jensen's Inequality.

```python
x = np.linspace(-10, 10, 1000)
y = x**2

f = 1 + 2*(x - 1)

_, ax = plt.subplots(figsize=(8,6))
ax.axhline(y=0, color='k')
ax.axvline(x=0, color='k')

ax.text(1.9, 2.2, ('$y = a + bx$'), color='xkcd:orange')
ax.text(1.22, 0.95, ('($x, g(x))$'))

ax.set_xticks([])
ax.set_yticks([])
ax.set_xticks([])
ax.set_yticks([])
ax.set_xlim([-2,4])
ax.set_ylim([-2,4])
ax.plot(x, y, 'xkcd:azure', lw=4, alpha=0.6)
ax.plot(x, f, 'xkcd:orange', linestyle='--', lw=4, alpha=0.6)

ax.title.set_text("line y tangent to g(x)")
```

Referring to the geometric figure above...

\begin{align}
  g(x) &\geq a + bx &\quad g(x) \text{ lies above the tangent } y = a + bx \\
  \Rightarrow g(X) &\geq a + bX &\quad \text{inequality applies to r.v. as well} \\
  \mathbb{E}\left( g(X) \right) &\geq \mathbb{E}\left( a + bX \right) &\quad \mathbb{E} \text{ as operator} \\\\
  \mathbb{E}\left( a + bX \right) &= a + b \, \mathbb{E}(X) \\
  &= a + b \, \mu &\quad \text{let } \mu = \mathbb{E}(X) \\
  &= g(\mu) \\
  &= g\left(\mathbb{E}(X)\right) \\\\
  \therefore \mathbb{E}\left( g(X) \right) &\geq g\left(\mathbb{E}(X)\right) &\quad \blacksquare
\end{align}

It is also possible to use Taylor expansions to prove Jensen's Inequality, but that is something for another time.

### (3) Markov's Inequality: $P(|X| \geq a) \leq \frac{\mathbb{E}|X|}{a}$ for any $a \gt 0$

The strength of Markov's Inequality is not in its accuracy of _approximation_, but rather in it simplicity and generality.

#### Proof

We will use the fundamental bridge to convert the probability of an event to the expected value of the indicator of said event.

\begin{align}
  I_{|X| \geq a} &\geq |X| &\quad \text{always true for } I_{|X| \geq a} = 0 \text{ or } I_{|X| \geq a} = 1 \\
  a \, \mathbb{E} \, I_{|X| \geq a} \geq \mathbb{E} \, |X| \\
  \mathbb{E} \, I_{|X| \geq a} \geq \frac{\mathbb{E} \, |X|}{a}
\end{align}

#### Intuition

Suppose that we have 100 people.

Q1: Is it possible that at least 95% of the people are younger than the mean age of the group?

A1: Yes, if one of those people is really, really old, then that person's age can raise the mean age above a large portion of group.

Q2: Is it possible that at least 50% of the people are older than *twice* the mean age?

A2: No, since that 50% will raise the mean age beyond that!

### (4) Chebyshev's Inequality: $P(|X-\mu| \leq \frac{Var(X)}{a^2}$ for $\mu = \mathbb{E}X, a \gt 0$

An alternate way of putting Chebyshev's Inequality is

\begin{align}
  P\left(|X-\mu| > c \, SD(X) \right) \leq \frac{1}{c^2} &\quad \text{for } c > 0 \\
\end{align}

Consider this: the probability that $X$ is more than 2 standard deviations from its mean is at most $\frac{1}{4}$. Recall the Normal distribution and the 68,95,98.7 percent rule? Compared to that rule, Chebyshev's Inequality provides a cruder but still correct relation.

#### Proof

\begin{align}
  P(|X-\mu| \ge a) &= P\left( (X-\mu)^2 \ge a^2 \right) \\
  &\leq \frac{\mathbb{E}(X-\mu)^2}{a^2} &\quad \text{by Markov's Inequality} \\
  &\leq \frac{Var(X)}{a^2} &\quad \blacksquare
\end{align}


text:
# MMA

- author: Dmytro Karabash
- categories: [MMA, diophantine equations]
- permalink: /temp/

Table of Contents

1  Problem 1. Does equation $x^{2}+y^{2}+z^{2}=2 x y z \tag{1}$ have any non-zero integer solutions?

1.1  Problem 1 Solution 1

1.1.1  Insight 1: consider mod 2

1.2  Problem 1 Solution 2

1.2.1  Insight 1: Vieta

2  Insight 2: Substitue

## Problem 1. Does equation $x^{2}+y^{2}+z^{2}=2 x y z \tag{1}$ have any non-zero integer solutions?

#collapse

### Problem 1 Solution 1

#### Insight 1: consider mod 2

Right hand side is even so there are two cases x,y,z are all even or one of them is even.

Case 1: If one of them is even then LHS(1) is 2 (mod 4) and RHS(1) is 0 (mod 4). Contradiction.

Case 2:If all of them are even divide the equation by the highest power of two that divides all of them, i.e. $$x=2^k x_0,y=2^k y_0,z=2^k z_0$$ and at least one of $x_0,y_0,z_0$ is odd, but then dividing equation by $2^{2k}$ we get
$$x_0^{2}+y_0^{2}+z_0^{2}=2^{k+1} x y z $$
and you can apply the argument from case 1 obtaining contradiction.

#collapse

> Note: when writting a proof you can skip extra circle of logic and go straight to case 2, hence the proof would just be:
> Let $k$ be highest power of 2 dividing $x$,$y$ and $z$, i.e $x=2^k x_0,y=2^k y_0,z=2^k z_0$  and at least one of $x_0,y_0,z_0$ is odd.  Dividing equation by $2^{2k}$ we obtain
$x_0^{2}+y_0^{2}+z_0^{2}=2^{k+1} x y z \tag{2}$ and since at least one of $x_0,y_0,z_0$ is odd, and RHS(2) is even exactly two of $x_0,y_0,z_0$ must be odd but then LHS(2) is 2 mod 4 and RHS(2) is 0 mod 4. Contradiction.

#collapse

### Problem 1 Solution 2

#### Insight 1: Vieta
\begin{align}
\mathbb 0  &= x^{2}-(2 y z) x +y^{2}+z^{2} \\
x_1+x_2& = 2yz\\
x_1x_2& = y^2+z^2
\end{align}

```python

```

#collapse

## Insight 2: Substitue
\begin{align}
x_1+x_2 & = 2yz  \\
x_1x_2  & = y^2 + z^2
\end{align}
which gives us
\begin{align}
x_2 & = 2yz - x_1 \\
x_2 & = \frac{y^{2}+z^{2}}{x_1} \\
\end{align}

```python

```


text:
---
layout: post
title:  "Monte-Carlo Integration made easy"
desc: "A practical Python explanation of Monte-Carlo integration"
long_desc: "A simple tutorial on Monte Carlo integration using Python. In this article, we'll cover how MC integration works, when to use MC integration and also how to use it efficiently via importance sampling, all in easy to read Python."
date: ###DATE
categories: [tutorial]
tags: [statistics]
loc: ###LOC
permalink: ###LINK
redirect_from: "/montecarlo"

math: true
---
Integrating a function is tricky. A lot of the time, the math is beyond us. Or beyond me, at the very least, and so I turn to my computer, placing the burden on its silent, silicon shoulders. And yet this isn't the end of it, because there are a host of ways to perform numerical integration.

Do we want the simple rectangle rule? The superior trapezoidal rule? Simpson's rule? Do we want to adaptively sample? How many dimensions is this in anyway - 1D, 2D, 3D... 100D?

Unfortunately, every algorithm listed above falls over at higher dimensionality, simply because most of them are based off a grid. If you have a 100 points in a grid, for a 1D integral, thats easy, 100 points. For a 2D grid, well now its 10 thousand cells. For a 3D grid, thats a million voxels. This is exponential scaling. This is bad news.

So instead we turn to the amazing algorithm of Monte-Carlo integration. Monte-Carlo here means its based off random numbers (yes, I'm glossing over a lot), and so we perform Monte-Carlo integration essentially by just taking the average of our function after evaluating it at some random points.

Let's just illustrate this with an example, starting with Simpson's rule. Let's integrate the super simple function:

$$ \int_0^{1.5\pi} \sin(x) \ dx $$

```python
# Remove
from base import *
plt.rcParams['axes.xmargin'] = 0
plt.rcParams['axes.ymargin'] = 0
plt.rcParams['axes.prop_cycle'] = (
    cycler(color=['#50d460', '#dbce1a', '#F77F00', '#FCBF49', '#EAE2B7']) +
    cycler(linestyle=['-', '--', ':', '-.', '-']))
```

```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import simps

def fn(x):
    return np.sin(x) + 1

xs = np.linspace(0, 1.5 * np.pi, 100)
ys = fn(xs)
area = simps(ys, x=xs)

plt.plot(xs, ys, label="Function")
plt.fill_between(xs, 0, ys, alpha=0.1)
plt.text(1, 0.75, f"Area from Simps is {area:0.3f}", fontsize=12)
plt.xlabel("x"), plt.legend();
plt.ylim(0, 2.2); ###REMOVE
```

Great, so how would we use Monte-Carlo integration to get another esimtate?

```python
plt.rcParams['axes.prop_cycle'] = (
    cycler(color=['#3ccfb4', '#dbce1a', '#F77F00', '#FCBF49', '#EAE2B7']) +
    cycler(linestyle=['-', '--', ':', '-.', '-']))
```

```python
width = 1.5 * np.pi - 0  # The width from 0 to 1.5pi
samples = np.random.uniform(low=0, high=width, size=1000000)
mc_area = fn(samples).mean() * width

plt.plot(xs, ys, label="Function")
plt.fill_between(xs, 0, ys, alpha=0.1)
plt.text(1, 0.75, f"Area from Simps is {area:0.3f}", fontsize=12)
plt.text(1, 0.5, f"Area from MC Integration is {mc_area:0.3f}", fontsize=12)
plt.xlabel("x"), plt.legend();
plt.ylim(0, 2.2); ###REMOVE
```

Alright, so let's dig into this a bit:

1. Why does this work?
2. What is `width`?
3. Why did I have to ask for a million samples!?!?

Lets start with \#1: **How does this work?**

Monte-Carlo integration is all about that [Law of Large Numbers](https://en.wikipedia.org/wiki/Law_of_large_numbers). To summarise the Wiki page, the LLN states that if you do an experiment over and over, the average of your experiment should converge to the expected value. Our experiment here is "sampling the function (uniformly)", so the LLN says if we keep sampling it, the average result should converge to the mean of the function. This should be intuitive - if you roll a fair 6-sided die a lot and take an average, you'd expect that you'd get around the same amount of each number, which would give you an average of 3.5.

Let's merge in **What is `width`** now. If we have the average of a function over some arbitrary $x$-domain, to get the area we need to factor in how big that $x$-domain is. We are essentially finding the area of a rectangle `width` wide, with an average height given by our samples.

If we want to be more formal about this, what we are doing is combining both our original function

$$ f(x) = \sin(x), $$

and the probability density function that describes how we draw our samples. In our case, this function is - in English - uniformly between $0$ and $1.5\pi$, and in mathematics:

$$p(x) =
    \begin{cases}
      \frac{1}{1.5\pi}, & \text{if}\ 0 < x < 1.5\pi \\
      0, & \text{otherwise}
    \end{cases}$$

The "width" comes in to our final result when you add the probability in to our equation:

$$\begin{align}\int_0^{1.5\pi} \sin(x)\ dx &= \int_0^{1.5\pi} \sin(x)\ \  \frac{1.5\pi}{1.5\pi} \ dx\\
&= \int_0^{1.5\pi}  1.5\pi\sin(x)   \ \ p(x) \ dx
\end{align}$$

Sorry for the math, but hopefully you can see that if we separate the equation so that we can get our sample function on the right, the `width` factor comes out naturally. Conceptually, it's easier to think of it using the rectangle analogy above, but that doesn't generalise too well.

Finally, **why did we need so many samples?**

The answer is that I wanted to make sure it agreed very well with the result from Simpsons' rule. Monte-Carlo integration has uncertainty, but you can quantify that:

$$ \text{error} = \frac{\sigma(x)}{\sqrt{N}}$$

where $\sigma$ is the standard deviation, $x$ is what we average (so really our samples times our width), and $N$ is the number of points. For us, the plot should really look like this:

```python
plt.rcParams['axes.prop_cycle'] = (
    cycler(color=['#1ca7e8', '#dbce1a', '#F77F00', '#FCBF49', '#EAE2B7']) +
    cycler(linestyle=['-', '--', ':', '-.', '-']))
```

```python
error = np.std(samples * width) / np.sqrt(samples.size)

plt.plot(xs, ys, label="Function")
plt.fill_between(xs, 0, ys, alpha=0.1)
plt.text(0.7, 0.75, f"Area from Simps is {area:0.3f}", fontsize=12)
plt.text(0.7, 0.5, f"Area from MC Integration is {mc_area:0.3f}±{error:0.3f}", fontsize=12)
plt.xlabel("x"), plt.legend();
plt.ylim(0, 2.2); ###REMOVE
```

Of course, Simpsons' rule has error too, let's not forget that!

## Importance sampling!

Importance sampling is the way that we can improve the accuracy of our estimates. It's conceptually simple - in the plot above, we could get better accuracy if we estimated the peak between 1 and 2 more throughly than if we estimated the area after 4 more thoroughly. So why are we uniformly sampling our distribution, when some areas are much more important??

For a super easy example, lets change the function. We now care about

$$\int_{-\infty}^{\infty} \frac{1 + x^2}{\sqrt{2\pi}} e^{-x^2/2}\ dx.$$

Uniformly sampling this would be crazy - how can we sample from $-\infty$ to $\infty$??? Instead, what we do is we look at the function and we separate it out. We say, "Hey, this looks like a polynomial times a normal distribution". Or more formally:

$$\int_{-\infty}^{\infty} (1+x^2) \ \ \mathcal{N}(0, 1)\ dx,$$

where $\mathcal{N}(0,1)$ is a normal distribution, centered at 0, with a width of 1. And just like before, we now have two parts - the first part to calculate, and the second part we can sample from. To do this, and then create a plot showing each sample, is simple:

```python
plt.rcParams['axes.prop_cycle'] = (
    cycler(color=['#32a840', '#dbce1a', '#F77F00', '#FCBF49', '#EAE2B7']) +
    cycler(linestyle=['-', '--', ':', '-.', '-']))
np.random.seed(2)
```

```python
# MC integration here
samples_2 = np.random.normal(size=1000)
fn_2 = 1 + samples_2 ** 2
area_2 = fn_2.mean()
error_2 = np.std(fn_2) / np.sqrt(fn_2.size)

# Simps integration here
def fn2(xs):
    return (1 + xs**2) * np.exp(-(xs**2)/2) / np.sqrt(2 * np.pi)
xs = np.linspace(-5, 5, 200)
ys = fn2(xs)
area_simps = simps(ys, x=xs)

# And of course, plotting here
plt.plot(xs, ys, label="Function", lw=3)
plt.fill_between(xs, 0, ys, alpha=0.1)
plt.text(-4.8, 0.5, f"MC Area={area_2:0.2f}±{error_2:0.2f}", fontsize=12)
plt.text(-4.8, 0.43, f"Simps Area={area_simps:0.2f}", fontsize=12)
plt.plot((samples_2, samples_2), ([0 for i in samples_2], [fn2(i) for i in samples_2]),
         c='#1c93e8', lw=0.2, ls='-', zorder=-1, alpha=0.5)
plt.xlabel("x"), plt.legend();
plt.ylim(0, 0.55); ###REMOVE
```

!!!main

Where each blue horiztonal line shows us one specific sample.
You can see that for us to get close to Simpons' rule we need far less samples, because we're sampling more efficiently.

When using importance sampling, note that you don't need to have a probability function you can sample with *perfectly* in your equation. You can put any PDF in (just like we did with the uniform distribution), and simply divide the original equation by that PDF. Imaging if we changed our function from above just a tiny bit:

$$\int_{-\infty}^{\infty} \frac{1 + x^2}{\sqrt{2\pi}} e^{-x^4/4}\ dx.$$

```python
def fn3(xs):
    return (1 + xs**2) * np.exp(-(xs**4)/4) / np.sqrt(2 * np.pi)

ys3 = fn3(xs)
plt.plot(xs, ys, label="Original")
plt.plot(xs, ys3, label="Modified to be $x^4/4$")
plt.legend(); plt.xlabel("x");
plt.ylim(0, 0.65); ###REMOVE
```

That's fine! We can still use that normal distribution from before, we just add it into the equation. Normally, your function will not be nice and analytic like the one we've tried to use, so we can state in general:

$$\begin{align} \text{integral} &= \int_{-\infty}^{\infty} f(x) dx\\
&= \frac{f(x)}{p(x)} \ p(x) \ dx\\
\end{align}$$

where $p(x)$ in our example will be the normal distribution.

```python
plt.rcParams['axes.prop_cycle'] = (
    cycler(color=['#e81c1c', '#e89a1c', '#F77F00', '#FCBF49', '#EAE2B7']) +
    cycler(linestyle=['-', '--', ':', '-.', '-']))
np.random.seed(2)
```

```python
from scipy.stats import norm

# MC integration here
x_samp = norm.rvs(size=2000)
p_of_x = norm.pdf(x_samp)
vals = fn3(x_samp) / p_of_x
area = vals.mean()
error = np.std(vals) / np.sqrt(vals.size)

# Simps integration here
xs = np.linspace(-5, 5, 200)
ys = fn3(xs)
area_simps = simps(ys, x=xs)

# And of course, plotting here
plt.plot(xs, ys, label="Function", lw=3)
plt.fill_between(xs, 0, ys, alpha=0.1)
plt.text(-4.8, 0.5, f"MC Area={area:0.2f}±{error:0.2f}", fontsize=12)
plt.text(-4.8, 0.43, f"Simps Area={area_simps:0.2f}", fontsize=12)
plt.plot((x_samp, x_samp), ([0 for i in x_samp], [fn3(i) for i in x_samp]),
         c='#e89a1c', lw=0.2, ls='-', zorder=-1, alpha=0.3)
plt.xlabel("x"), plt.legend();
plt.ylim(0, 0.65); ###REMOVE
```

So hopefully you can see how this would be useful. To summarise, the general process for Monte-Carlo integration is:

1. Have your function to integrate. 1D, 2D, 3D, doesn't matter.
2. Pick a distribution that looks *something* like your function, or just use a Uniform distribution
3. Sample points from that distribution
4. Get the function at those points, and divide by $p(x)$.
5. Take the mean for the estimate, and the standard deviation / root(N) for the error.
6. Celebrate

Finally, obviously I've kept the examples here to 1D for simplicity, but I really should stress that MC integration shines in higher dimensions. If you have a 10 dimensional function that looks roughly Gaussian (like a normal), you can sample from a 10 dimensional normal, and apply all the same steps above, nothing at all changes.

max_stars_repo_name: Samreay/samreay

max_issues_repo_name: Samreay/samreay

max_forks_repo_name: Samreay/samreay


text:
# 3M1 Examples Paper crib (linear algebra)

Garth N. Wells (gnw20@cam.ac.uk)

(C) 2016-2020

```python
import numpy as np
from scipy import linalg
```

## Question 1

The usual rules for real matrices extend without modification to complex valued matrices.

a. $\det(\boldsymbol{A}) = (4 + 4i)(4 + i) - (2-i)(-3 + 2i) = 16 + 13 i$

b. $\det(\boldsymbol{A}^{H}) = (4 - 4i)(4 - i) - (-3 - 2i)(2 + i) = 16 - 13i$

c. Inverse:

   $$
   \boldsymbol{A}^{-1} = \frac{1}{16 + 13 i }
   \begin{bmatrix}
      4 + i   & -2 + i  \\
      3 - 2i  & 4 + 4i
   \end{bmatrix}
   $$

**Extra:** [SymPy](https://www.sympy.org/) can be used for symbolic computations:

```python
import sympy
sympy.interactive.printing.init_printing(use_unicode=False, wrap_line=False)
A = sympy.Matrix([[4 + sympy.I*4 , 2- sympy.I], [-3 + sympy.I*2 , 4 + sympy.I]])
A
```

```python
det = A.det()
sympy.simplify(det)
```

```python
detH = A.H.det()
sympy.simplify(detH)
```

```python
Ainv = A.inv()
sympy.simplify(Ainv)
```

```python
sympy.simplify(A.multiply(A.inv()))
```

## Question 2

The diagonal entries must be real since $A_{ii} = \bar{A}_{ii}$.

## Question 3

\begin{align}
\det \left(\boldsymbol{Q} - \lambda \boldsymbol{I} \right) &=
\det\left(
\begin{bmatrix}
\cos \theta - \lambda & - \sin \theta \\
\sin \theta           & \cos \theta - \lambda
\end{bmatrix} \right)
\\
&= (\cos \theta - \lambda)^{2} + \sin^{2}\theta \\
&= \lambda^{2} - (2 \cos \theta) \lambda  + 1 \\
&= 0
\end{align}

Computing roots,

\begin{align}
\lambda &= \cos \theta \pm \sqrt{\cos^{2}\theta - 1} \\
        &= \cos \theta \pm i \sin\theta
\end{align}

## Question 4

a. If $\boldsymbol{x}$ is an eigenvector and $\lambda$ is an eigenvalue,

   $$
   \boldsymbol{M} \boldsymbol{x} = \lambda \boldsymbol{x}
   $$

   Premultiplying by $\boldsymbol{x}^{H}$,

   $$
   \boldsymbol{x}^{H} \boldsymbol{M} \boldsymbol{x} = \lambda \boldsymbol{x}^{H} \boldsymbol{x}
   $$

   Taking the complex conjugate of both sides and noting that $\boldsymbol{x}^{H} \boldsymbol{x}$ is real,

   $$
   \left(\boldsymbol{x}^{H} \boldsymbol{M} \boldsymbol{x}\right)^{H} =
   \boldsymbol{x}^{H} \boldsymbol{M}^{H} \boldsymbol{x} =
   \bar{\lambda} \boldsymbol{x}^{H} \boldsymbol{x}
   $$

   Since $\boldsymbol{M}$ is Hermitian ($\boldsymbol{M} =\boldsymbol{M}^{H}$),
   comparing the above two equations they can hold ony if $\lambda = \bar{\lambda}$,
   which is true only if $\lambda$ is *real*.

b. Consider two eigenpairs $(\lambda_{1}, \boldsymbol{x}_{1})$ and $( \lambda_{2}, \boldsymbol{x}_{2})$, where $\lambda_{1} \ne \lambda_{2}$. We have:

   \begin{align}
   \boldsymbol{M} \boldsymbol{x}_{1} &= \lambda_{1} \boldsymbol{x}_{1} \\
   \boldsymbol{M} \boldsymbol{x}_{2} &= \lambda_{2} \boldsymbol{x}_{2}
   \end{align}

   Multiplying by $\boldsymbol{x}_{2}$ and $\boldsymbol{x}_{1}$, respectively,

   \begin{align}
   \boldsymbol{x}_{2}^{H} \boldsymbol{M} \boldsymbol{x}_{1} &= \lambda_{1} \boldsymbol{x}_{2}^{H}
   \boldsymbol{x}_{1} \\
   \boldsymbol{x}_{1}^{H} \boldsymbol{M} \boldsymbol{x}_{2} &= \lambda_{2} \boldsymbol{x}_{1}^{H}
   \boldsymbol{x}_{2}
   \end{align}

   Taking the complex conjugate transpose of the first equation and exploiting that $\boldsymbol{M} =
   \boldsymbol{M}^{H}$,

   \begin{equation}
   \left(\boldsymbol{x}_{2}^{H} \boldsymbol{M} \boldsymbol{x}_{1}\right)^{H} =  \boldsymbol{x}_{1}^{H}
   \boldsymbol{M} \boldsymbol{x}_{2} = \lambda_{1} \boldsymbol{x}_{1}^{H} \boldsymbol{x}_{2}
   \end{equation}

   Comparing this to $\boldsymbol{x}_{1}^{H}\boldsymbol{M} \boldsymbol{x}_{2} = \lambda_{2}
   \boldsymbol{x}_{1}^{H}\boldsymbol{x}_{2}$, since $\lambda_{1} \ne
   \lambda_{2}$ both equations can hold only if $\boldsymbol{x}_{1}^{H} \boldsymbol{x}_{2} = 0$, i.e. the   eigenvectors are *orthogonal*.

## Question 5

To show that $\boldsymbol{A}^{H} \boldsymbol{A}$ is positive semi-definite:
\begin{equation}
\boldsymbol{x}^{H} \boldsymbol{A}^{H} \boldsymbol{A} \boldsymbol{x}  =
(\boldsymbol{A}\boldsymbol{x})^{H} (\boldsymbol{A} \boldsymbol{x})
\ge 0
\end{equation}
If $(\lambda, \boldsymbol{x})$ is an eigenpair of $\boldsymbol{A}^{H} \boldsymbol{A}$,

\begin{equation}
  \boldsymbol{A}^{H} \boldsymbol{A} \boldsymbol{x}
  = \lambda \boldsymbol{x}
\end{equation}

Premultiplying both sides by $\boldsymbol{x}^{H}$:

\begin{equation}
  \boldsymbol{x}^{H}  \boldsymbol{A}^{H} \boldsymbol{A} \boldsymbol{x}
  = \lambda \boldsymbol{x}^{H} \boldsymbol{x}
\end{equation}

Since $\boldsymbol{x}^{H} \boldsymbol{A}^{H} \boldsymbol{A} \boldsymbol{x} \ge 0$ and $\boldsymbol{x}^{H} \boldsymbol{x} \ge 0$, all eigenvalues $\lambda$ must be greater than or equal to zero.

## Question 6

a. Eigenvalues of $\boldsymbol{A}$ satisfy $\det(\boldsymbol{A} - \lambda\boldsymbol{I}) = 0$.
   For the transpose, we have
   \begin{equation}
     \det\left(\boldsymbol{A}^{T} - \lambda\boldsymbol{I}\right) =
     \det\left(\left(\boldsymbol{A} - \lambda \boldsymbol{I}\right)^{T}\right) =
     \det(\boldsymbol{A} - \lambda \boldsymbol{I}),
   \end{equation}
   since the $\det\boldsymbol{A} = \det \boldsymbol{A}^{T}$. Hence eigenvalues of $\det\boldsymbol{A}$ and
   $\det \boldsymbol{A}^{T}$ are the same.

b. Noting that $\det\left(\boldsymbol{A}^{H}\right) = \det\left(\overline{\boldsymbol{A}}^{T}\right) = \overline{\det\left(\boldsymbol{A}^{T}\right)} = \overline{\det\left(\boldsymbol{A}\right)}$,
   \begin{equation}
   \det\left(\boldsymbol{A}^{H} - \lambda\boldsymbol{I}\right) =
   \det\left(\left(\boldsymbol{A} - \bar{\lambda} \boldsymbol{I}\right)^{H}\right) =
   \det\left( \overline{(\boldsymbol{A} - \bar{\lambda}\boldsymbol{I})^{T}}\right) =
   \overline{\det(\boldsymbol{A} - \bar{\lambda} \boldsymbol{I}}).
   \end{equation}

c. If $\lambda$ and $\boldsymbol{x}$ are an eigenvalue and eigenvector, respectively, of
   $\boldsymbol{A}\boldsymbol{B}$, then
   \begin{equation}
    \boldsymbol{A}\boldsymbol{B} \boldsymbol{x} = \lambda \boldsymbol{x}
   \end{equation}
   Multiplying $\boldsymbol{x}$ by $\boldsymbol{B}$ on both sides,
   \begin{equation}
     (\boldsymbol{B}\boldsymbol{A}) (\boldsymbol{B} \boldsymbol{x})
     = \lambda (\boldsymbol{B}\boldsymbol{x})
   \end{equation}
   We see that $\lambda$ is an eigenvalue of $\boldsymbol{B}\boldsymbol{A}$, and the eigenvector is now
   $\boldsymbol{B} \boldsymbol{x}$.

## Question 7

Follows directly from the definition of the norms:

\begin{equation}
  \max_{i=1}^{n} |x_{i}| \le \sum_{i=1}^{n} |x_{i}| \le n \max_{i=1}^{n} |x_{i}|.
\end{equation}

## Question 8

Recall that from the defintion of a matrix operator norm it follows that $\| \boldsymbol{A}\boldsymbol{x} \| \le \| \boldsymbol{A} \| \| \boldsymbol{x} \|$

a. From the definition of the matrix norm:

   \begin{equation}
      \| \boldsymbol{A} \boldsymbol{B} \boldsymbol{x} \|
      \le \| \boldsymbol{A} \| \| \boldsymbol{B} \boldsymbol{x}\|
      \le \| \boldsymbol{A} \| \| \boldsymbol{B}\| \| \boldsymbol{x}\|
   \end{equation}

   for all $\boldsymbol{x}$. Re-arranging,

   \begin{equation}
      \frac{\| \boldsymbol{A} \boldsymbol{B} \boldsymbol{x} \|}{\|\boldsymbol{x}\|}
      \le \| \boldsymbol{A} \| \| \boldsymbol{B}\|
   \end{equation}

   for all $\boldsymbol{x} \ne \boldsymbol{0}$. From the definition of the norm, $\|
   \boldsymbol{A} \boldsymbol{B} \|$ is the largest possible value of $\| \boldsymbol{A}
   \boldsymbol{B} \boldsymbol{x} \| / \|\boldsymbol{x}\|$ (over all
   $\boldsymbol{x}$, exluding the zero vector), hence

   \begin{equation}
      \| \boldsymbol{A} \boldsymbol{B} \|  \le \| \boldsymbol{A} \| \| \boldsymbol{B}\|
   \end{equation}

b. Consider
   \begin{equation}
      \| (\boldsymbol{A} + \boldsymbol{B}) \boldsymbol{x} \|
      = \| \boldsymbol{A}\boldsymbol{x} + \boldsymbol{B} \boldsymbol{x} \|
   \end{equation}

   From the triagle inequality for vectors,

   \begin{align}
      \| \boldsymbol{A}\boldsymbol{x} + \boldsymbol{B} \boldsymbol{x} \|
      &\le
      \| \boldsymbol{A}\boldsymbol{x} \| + \| \boldsymbol{B} \boldsymbol{x} \|  \\
      &\le
      \| \boldsymbol{A} \| \| \boldsymbol{x} \| + \| \boldsymbol{B} \| |\boldsymbol{x} \|
   \end{align}

   Re-arranging

   \begin{equation}
      \frac{\| (\boldsymbol{A} + \boldsymbol{B}) \boldsymbol{x} \|}{\| \boldsymbol{x} \|}
      \le
      \| \boldsymbol{A} \| + \| \boldsymbol{B} \|
   \end{equation}

   Using same argument as from part a), we get $\| \boldsymbol{A} + \boldsymbol{B} \| \le
   \| \boldsymbol{A} \| + \| \boldsymbol{B} \|$.

## Question 9

Recall that the smallest possible value of $C$ for which $\|{\boldsymbol{A}\boldsymbol{x}}\| \le C \| \boldsymbol{x}\|$ holds is the operator norm $\|\boldsymbol{A}\|$. The task is to find $C$ for the different norms. It simplifies the proofs if we consider all vectors for which $\|\boldsymbol{x}\| = 1$ (this is possible since $\|\alpha \boldsymbol{x}\| = |\alpha| \|\boldsymbol{x}\|$), in which case the task is to find the smallest possible $C$ such that $\|\boldsymbol{A}\boldsymbol{x}\| \le C$.

What we need to do is (i) prove the inequality, and then (ii) find an equality to show that is is a weak inequality.

a. For the 1-norm, denoting the $j$th column of $\boldsymbol{A}$ by
   $\boldsymbol{a}_{j}$, and for a vector $\|\boldsymbol{x}\|_{1} = 1$:

   \begin{align*}
        \| \boldsymbol{A} \boldsymbol{x}\|_{1}
        = \|\sum_{j=1}^{n} \boldsymbol{a}_{j} x_{j}\|_{1}
        \le
        \sum_{j=1}^{n} \|\boldsymbol{a}_{j}\|_{1} |x_{j}|
        \le \max_{j=1}^{n} \|\boldsymbol{a}_{j}\|_{1}
   \end{align*}

   This is the maximum column sum.

   The term on the right could possibly be larger than $\| \boldsymbol{A} \|_{1}$, whereas the norm is the smallest
   possible value for the RHS that still satisfies the inequalities.  If we can show a case for which equality is
   reached, $\|\boldsymbol{A} \boldsymbol{x}\|_{1} = \max_{j=1}^{n} \|\boldsymbol{a}_{j}\|_{1}$, we have the norm.
   For the vector the $\boldsymbol{e}_{j}$ with $e_{j} = 1$, where $j$ is the column with the greatest 1-norm,
   and $e_{i \ne j} = 0$, we have
   \begin{equation}
        \|\boldsymbol{A} \boldsymbol{e}_{j}\|_{1} = \max_{j=1}^{n} \|\boldsymbol{a}_{j}\|_{1}
   \end{equation}
   Therefore,
   \begin{equation}
      \|\boldsymbol{A}\|_{1} = \max_{j=1}^{n} \|\boldsymbol{a}_{j}\|_{1}
   \end{equation}

b. For a vector $\|\boldsymbol{x}\|_{\infty} = 1$:

   \begin{align*}
        \|\boldsymbol{A} \boldsymbol{x}\|_{\infty}
          = \max_{i=1}^{m} |\sum_{j=1}^{n} a_{ij} x_{j}|
          &\le \max_{i=1}^{m} \sum_{j=1}^{n} |a_{ij}| |x_{j}|
          \\
          &\le \max_{i=1}^{m} \sum_{j=1}^{n} |a_{ij}|
   \end{align*}

   The RHS is the maximum row sum.

   As before, we need to find a case with equality.  If the row
   with the maximum sum is row $k$, we choose a vector $\boldsymbol{x}$
   where $x_{i} = \pm 1$ such that the sign of $x_{i}$ is the
   same as the sign of the entry $a_{ki}$. We then have

   \begin{equation}
      \|\boldsymbol{A}\|_{\infty}
      = \max_{i=1}^{m} \sum_{j=1}^{n} |a_{ij}|
   \end{equation}

## Question 10

Recall that the $l_{2}$ norm is the square root of the maximum eigenvalue of $\boldsymbol{A}^{T}\boldsymbol{A}$ (matrix is real in this question). We have
$$
\boldsymbol{A}^{T}\boldsymbol{A} =
\begin{bmatrix}
26 & 5 \\ 5 & 1
\end{bmatrix}
$$
For this matrix, $\lambda = (27 \pm \sqrt{27^{2} -4})/2$. Hence $\|\boldsymbol{A}\|_{2} \approx 5.1926$.

Checking for a vector of length $1$:
$$
\boldsymbol{A}
\begin{bmatrix} \cos \theta \\ \sin \theta   \end{bmatrix} =
\begin{bmatrix} \cos \theta \\ 5\cos \theta + \sin \theta  \end{bmatrix}
$$
Comptuing the norm,
$$
F(\theta) = \|\boldsymbol{A}\boldsymbol{x}\|_{2}^{2} = 25 \cos^{2} \theta + 10 \cos\theta\sin\theta  + 1
$$

To find the extreme points of the function, we differentiate $F$ with respect to $\theta$  and set the derivative equal to zero:
$$
\begin{align*}
\frac{d F}{d \theta} &= -50 \cos\theta \sin\theta + 10 (\cos^{2} \theta - \sin^{2} \theta) \\
&= -25 \sin 2\theta + 10 \cos 2\theta \\
&= 0
\end{align*}
$$
The norm is maximised when $\theta = \arctan(2/5) /2$. Evaluating:

```python
θ = np.arctan(2/5)/2
print("θ =", θ)

F = 25*np.cos(θ)*np.cos(θ) + 10*np.cos(θ)*np.sin(θ) + 1
print("F =", np.sqrt(F))
```

    θ = 0.19025318855618245
    F = 5.192582403567252

## Question 11 (condition number)

Recall that $\kappa(\boldsymbol{A}) = \| \boldsymbol{A} \| \| \boldsymbol{A}^{-1} \|$.

We have $\| \boldsymbol{A}\|_{1} = 3$ and $\| \boldsymbol{A}\|_{\infty} \approx 2$.
Note that
$$
\boldsymbol{A}^{-1} \approx -\frac{1}{2}
\begin{bmatrix}
1 & -2 \\
-1 & 10^{-4}
\end{bmatrix}
$$
so $\| \boldsymbol{A}^{-1} \|_{1} \approx 1$ and $\| \boldsymbol{A}^{-1} \|_{\infty} = 3/2$.
Therefore $\kappa_{1}(\boldsymbol{A}) \approx 3$ and $\kappa_{\infty}(\boldsymbol{A}) \approx 3$.

Recall that $\kappa_{2}(\boldsymbol{A}) = \sqrt{\lambda_{\max}(\boldsymbol{A}^{T}\boldsymbol{A}})/\sqrt{\lambda_{\min}(\boldsymbol{A}^{T}\boldsymbol{A}})$. As an approximation we ignore the $(1, 1)$ entry:
$$
  \boldsymbol{A}^{T}\boldsymbol{A} \approx
  \begin{bmatrix}
     1 & 1 \\ 1 & 5
  \end{bmatrix}
$$
Computing the eigenvalues, $\lambda \approx 3 \pm \sqrt{5}$, therefore $\kappa_{2} \approx \sqrt{3 + \sqrt{5}}/\sqrt{3 - \sqrt{5}} \approx 2.618$.

The matrix is very well conditioned, however LU factorisation may require pivoting. This is an issue with LU factorisation rather than a pathological problem with the matrix.

## Question 12 (least squares)

A projection matrix $\boldsymbol{P}$  has the property $\boldsymbol{P} = \boldsymbol{P}\boldsymbol{P}$, and $\boldsymbol{P}^{H} = \boldsymbol{P}$.

a. The solution to the least squares problem is $\hat{\boldsymbol{x}} = \boldsymbol{A}
   (\boldsymbol{A}^{H} \boldsymbol{A})^{-1} \boldsymbol{A}^{H} \boldsymbol{b}$.
   Therefore $\boldsymbol{r} = \boldsymbol{A} \hat{\boldsymbol{x}} - \boldsymbol{b}
   = \boldsymbol{A}(\boldsymbol{A}^{H} \boldsymbol{A})^{-1} \boldsymbol{A}^{H}
   \boldsymbol{b} - \boldsymbol{b}$. Insert this expression for $\boldsymbol{r}$ into the expression in the
   question and re-arrange. to show the result.

   Vectors $\boldsymbol{A}\boldsymbol{z}$ lie in the *column space* of $\boldsymbol{A}$,
   hence the expression says that the least-squares residual is *orthogonal* to the column
   space of $\boldsymbol{A}$.

b.
   $$
     \boldsymbol{P}\boldsymbol{P}
        = \boldsymbol{A} (\boldsymbol{A}^{H}\boldsymbol{A})^{-1}\boldsymbol{A}^{H}
    \boldsymbol{A}(\boldsymbol{A}^{H}\boldsymbol{A})^{-1}\boldsymbol{A}^{H}
     = \boldsymbol{A} (\boldsymbol{A}^{H}\boldsymbol{A})^{-1}\boldsymbol{A}^{H}
    = \boldsymbol{P}
   $$
   and
   $$
        \boldsymbol{P}^{H}
        = \boldsymbol{A}(\boldsymbol{A}^{H}\boldsymbol{A})^{-H}\boldsymbol{A}^{H}
        = \boldsymbol{A}(\boldsymbol{A}^{H}\boldsymbol{A})^{-1}\boldsymbol{A}^{H}
   $$
   by $\boldsymbol{A}^{H}\boldsymbol{A}$ being Hermitian

c. We can phrase a least squares problem as
   $$
        \boldsymbol{A} \hat{\boldsymbol{x}}
        = \boldsymbol{A}(\boldsymbol{A}^{H}\boldsymbol{A})^{-1}\boldsymbol{A}^{H}
        \boldsymbol{b}
        = \boldsymbol{P} \boldsymbol{b}
   $$
   which says that $\boldsymbol{P}$ projects $\boldsymbol{b}$ into the column space of $\boldsymbol{A}$. If
   $\boldsymbol{b}^{\prime} = \boldsymbol{P}\boldsymbol{b}$ is in the column space of $\boldsymbol{A}$, then
   $\boldsymbol{A} \hat{\boldsymbol{x}} = \boldsymbol{b}^{\prime}$ has a solution.

## Question 13 (pseudo inverse)

a.

   - Firstly, if the $m \times n$ matrix $\boldsymbol{A}$ has linearly independent
     rows, then the rank of $\boldsymbol{A}$ is $m$ and the column space of
     $\boldsymbol{A}$ spans $\mathbb{C}^{m}$, and the nullspace space
     of $\boldsymbol{A}^{H}$ contains the zero vector only.

   - Now, consider the nullspace of $\boldsymbol{A}\boldsymbol{A}^{H}$:

     \begin{equation}
      \boldsymbol{A}\boldsymbol{A}^{H} \boldsymbol{x} = \boldsymbol{0} \ \rightarrow
      \ \boldsymbol{x}^{H}\boldsymbol{A} \boldsymbol{A}^{H} \boldsymbol{x} = 0 \ \rightarrow
      \ (\boldsymbol{A}^{H} \boldsymbol{x})^{H}\boldsymbol{A}^{H} \boldsymbol{x} = 0
     \end{equation}

     The above holds only if $\boldsymbol{A}^{H} \boldsymbol{x} = \boldsymbol{0}$, which
     says that $\boldsymbol{x}$ must come from the nullspace of
     $\boldsymbol{A}^{H}$.  We have already determined that the nullspace of
     $\boldsymbol{A}^{H}$ contains only the zero vector, therefore
     $\boldsymbol{A}\boldsymbol{A}^{H}$ is full rank (the nullspace of
     $\boldsymbol{A}\boldsymbol{A}^{H}$ contains the zero vector only) and can be inverted.

b. Since $\boldsymbol{A}\boldsymbol{A}^{H}$ is square and full rank, it can be inverted,
   $$
     \boldsymbol{A} \boldsymbol{A}^{+} = \boldsymbol{A} \boldsymbol{A}^{H}
     (\boldsymbol{A}\boldsymbol{A}^{H})^{-1} = \boldsymbol{I}.
   $$
   hence $\boldsymbol{A}^{+}$ is a right-inverse of $\boldsymbol{A}$.

## Question 14 (stationary iterative methods)

Define matrix $\boldsymbol{A}$:

```python
A = np.array([[2, -1], [-1, 2]])
print(A)
```

    [[ 2 -1]
     [-1  2]]

We split $\boldsymbol{A}$ such that $\boldsymbol{A} = \boldsymbol{N} - \boldsymbol{P}$. A method will converge if the largest absolute eigenvalue of $\boldsymbol{N}^{-1}\boldsymbol{P}$ is less the one.

For the Richardson method $\boldsymbol{N} = \boldsymbol{I}$:

```python
# Richardson
N = np.identity(2)
P = N - A
M = np.linalg.inv(N).dot(P)
print(np.linalg.eigvals(M))
```

    [ 0. -2.]

Largest eigenvalue (absolute value) is greater than 1, therefore method will not converge.

For the Jacobi method, $\boldsymbol{N} = \text{diag}(\boldsymbol{A})$:

```python
# Jacobi
N = np.diag(np.diag(A))
P = N - A
M = np.linalg.inv(N).dot(P)
print(np.linalg.eigvals(M))
```

    [ 0.5 -0.5]

Largest eigenvalue (absolute value) is less than 1, therefore method will converge.

For Gauss-Seidel, $\boldsymbol{N}$ is the lower triangular part of $\boldsymbol{A}$:

```python
# Gauss-Seidel
N = np.tril(A)
P = N - A
M = np.linalg.inv(N).dot(P)
print(np.linalg.eigvals(M))
```

    [0.   0.25]

Gauss-Seidel will converge because largest eigenvalue is less than one.

## Question 15 (SVD)

\begin{align}
  \boldsymbol{A}^{-1} &= \left(\boldsymbol{U} \boldsymbol{\Sigma} \boldsymbol{V}^{H}\right)^{-1} \\
  &=\boldsymbol{V}^{-H} \boldsymbol{\Sigma}^{-1} \boldsymbol{U}^{-1} \\
  &=\boldsymbol{V} \boldsymbol{\Sigma}^{-1} \boldsymbol{U}^{H}
\end{align}
Non-singular matrix cannot have any zero singular values. In fact, smallest singular values is a measure of the 'distance' to a singular matrix.

## Question 16 (SVD)

Define matrix $\boldsymbol{A}$:

```python
A = np.array([[1, 3], [2, 2], [3, 1]])
print(A)
```

    [[1 3]
     [2 2]
     [3 1]]

Compute the reduced SVD (recall that NumPy uses $\boldsymbol{A} = \boldsymbol{U} \boldsymbol{\Sigma} \boldsymbol{V}$ rather than  $\boldsymbol{A} = \boldsymbol{U} \boldsymbol{\Sigma} \boldsymbol{V}^{T}$):

```python
U, s, V = np.linalg.svd(A, full_matrices=False)
print(s)
print(U)
print(V)
```

    [4.89897949 2.        ]
    [[-5.77350269e-01  7.07106781e-01]
     [-5.77350269e-01  8.98662938e-17]
     [-5.77350269e-01 -7.07106781e-01]]
    [[-0.70710678 -0.70710678]
     [-0.70710678  0.70710678]]

Compute pseudoinverse by creating $\boldsymbol{\Sigma}^{+} =  \boldsymbol{\Sigma}_{1}^{-1}$

```python
# Pseudoinverse
Sigma_p = np.diag(1.0/s)
Ap = (V.T).dot(Sigma_p.dot(U.T))
print(Ap)
```

    [[-0.16666667  0.08333333  0.33333333]
     [ 0.33333333  0.08333333 -0.16666667]]

Compute $\boldsymbol{A}^{+}\boldsymbol{A}$:

```python
# Check that A^{+}A = I
print(Ap.dot(A))
```

    [[ 1.00000000e+00 -3.88578059e-16]
     [ 2.49800181e-16  1.00000000e+00]]

which is the identity. Compute now $\boldsymbol{A}\boldsymbol{A}^{+}$

```python
# Check that AA^{+} \ne I
print(A.dot(Ap))
```

    [[ 0.83333333  0.33333333 -0.16666667]
     [ 0.33333333  0.33333333  0.33333333]
     [-0.16666667  0.33333333  0.83333333]]

Which is clearly not the identity.

(e) Recall that from $\boldsymbol{A}^{T} \boldsymbol{A} \hat{\boldsymbol{x}} = \boldsymbol{A}^{T} \boldsymbol{b}$ we have
\begin{align}
\hat{\boldsymbol{x}} &= (\boldsymbol{A}^{T} \boldsymbol{A})^{-1}\boldsymbol{A}^{T} \boldsymbol{b} \\
&= \boldsymbol{A}^{+} \boldsymbol{b}
\end{align}
Multiplying both sides by $\boldsymbol{A}$,
$$
\boldsymbol{A}\hat{\boldsymbol{x}} = \underbrace{\boldsymbol{A} \boldsymbol{A}^{+}}_{\boldsymbol{P}} \boldsymbol{b}
$$
where $\boldsymbol{P}$ is the projection matrix from an earlier questions. Recall that $\boldsymbol{P}$ projects a vector into the column space of $\boldsymbol{A}$.
Since $\boldsymbol{P}$ is a projector, it does nothing if $\boldsymbol{b}$ is already in column space. Therefore any vector $\boldsymbol{b}$ in column space of $\boldsymbol{A}$ is a solution.

## Question 17 (pseudo inverse)

Define matrix $\boldsymbol{A}$:

```python
A = np.array([[0, 3, 0], [2, 0, 0]])
print(A)
```

    [[0 3 0]
     [2 0 0]]

Compute SVD:

```python
U, s, V = np.linalg.svd(A, full_matrices=False)
print(U)
print(s)
print(V)
```

    [[1. 0.]
     [0. 1.]]
    [3. 2.]
    [[0. 1. 0.]
     [1. 0. 0.]]

Compute $\boldsymbol{\Sigma}^{+} =  \boldsymbol{\Sigma}_{1}^{-1}$, and then $\boldsymbol{A}^{+}$

```python
# \Sigma^{+}
Sigma_p = np.diag(1.0/s)

# A^{+}
Ap = (V.T).dot(Sigma_p).dot(U)
print(Ap)
```

    [[0.         0.5       ]
     [0.33333333 0.        ]
     [0.         0.        ]]

Compute $\boldsymbol{A}^{+}\boldsymbol{A}$

```python
ApA = Ap.dot(A)
print(ApA)
```

    [[1. 0. 0.]
     [0. 1. 0.]
     [0. 0. 0.]]

For this matrix, any $\boldsymbol{x} = [x_{1} \ \ x_{2} \ \ 0]$ (which is from the row space of $\boldsymbol{A}$) satisfies $\boldsymbol{A}^{+} \boldsymbol{A} \boldsymbol{x} = \boldsymbol{x}$.

## Question 18 (rank deficient least squares)

Define matrix and RHS vector

```python
A = np.array([[1, 0, 0], [1, 0, 0], [1, 1, 1]])
print(A)

b = np.array([0, 2, 2])
print(b)
```

    [[1 0 0]
     [1 0 0]
     [1 1 1]]
    [0 2 2]

Compute the SVD and print singular values:

```python
# Compute SVD
U, s, V = np.linalg.svd(A)
print(s)
```

    [2. 1. 0.]

There is one zero singular value, so we need to 'trim' the last column from $\boldsymbol{U}$ and the last row from $\boldsymbol{U}$, and compute $\boldsymbol{\sigma}^{+}$

```python
# Create view of U with last columns
U1 = U[:, :2]
print(U1)
V1 = V[:2,:]
print(V1)

# Create Sigma^{+}
S1 = np.diag(1.0/s[:-1])
print(S1)
```

    [[-0.40824829 -0.57735027]
     [-0.40824829 -0.57735027]
     [-0.81649658  0.57735027]]
    [[-0.81649658 -0.40824829 -0.40824829]
     [-0.57735027  0.57735027  0.57735027]]
    [[0.5 0. ]
     [0.  1. ]]

Solve problem $\boldsymbol{x} = \boldsymbol{V}_{1} \boldsymbol{\Sigma}_{1}^{-1} \boldsymbol{U}_{1}^{T}$

```python
x = np.transpose(V1).dot(S1.dot(U1.T).dot(b))
print(x)
```

    [1.  0.5 0.5]

max_stars_repo_path: 3M1 Examples Paper 1 Crib.ipynb

max_issues_repo_path: 3M1 Examples Paper 1 Crib.ipynb

max_forks_repo_path: 3M1 Examples Paper 1 Crib.ipynb


text:
# Physics 420/580 Midterm Exam
## October 19, 2017 1pm-2pm

Do the following problems.  Use the Jupyter notebook, inserting your code and any textual answers/explanations in cells between the questions. (Feel free to add additional cells!)  Marks will be given based on how clearly you demonstrate your understanding.

There are no restrictions on downloading from the internet, eclass,  or the use of books, notes, or any other widely available computing resources.  However, **you are not allowed** to communicate with each other or collaborate in any way and uploading to the internet or sending or receiving direct communications is not appropriate.

When you are finished, upload the jupyter notebook to eclass. Eclass times out after 2:05 so make sure that you upload things before then.  Also be careful to save the notebook periodically and that you upload your final exam file.

```python

```

## Graphics

Plot the two curves:
$$\begin{align} y&=4x^3-3x-2\\
x&=\sin(\frac{y^4}{4}-2y^2+2)
\end{align}
$$

How many intersections are there?  Read the x,y value corresponding to the intersections from the plots.

```python
import numpy as np
import matplotlib.pyplot as plt
x=np.linspace(-1.2,1.2,1000)
y=np.linspace(-3.5,1,1000)
plt.plot(x,4*x**3-3*x-2)
plt.plot(np.sin(y**4/4-2*y**2+2),y)
plt.grid()
plt.show()

```

There are six intersections (-.9,-2.4), (-.75,-1.3), (-.3,-1.2), (.25,-2.6), (.35,-2.8), (.65,-2.8)(.7,-2.6)

# Reconstruction

A instantaneous flash of light occurs in a large tank of water at time $t_0$ and at position $\vec{x_0}$

The group velocity of light in water is about $2.2 \times 10^8$ m/s.  Four sensors detect the light flash, and report a measurement of the time at which light strikes the sensor.  The locations of the sensors and the time which each sensor was hit is recorded in the table below:

|Sensor #|	x|	y|	z|		Time|
|--------|------------|------------|------------|------------|
|	|[m]|	[m]|	[m]|		[s]|
|1	|0	|0|	10	|	6.3859E-08|
|2	|8.66025404|	0|	-5|		1.1032E-07|
|3	|-4.33012702|	7.5|	-5|		7.9394E-08|
|4	|-4.33012702|	-7.5|	-5|		1.0759E-07|

Calculate the initial time $t_0$ and the location of the flash.

Here is a sketch of the general geometry.  Note that the gray region is not different from the rest of the water- it is just to show the tetrahedral arrangement of the light sensors.  The flash, shown as the yellow star is in an arbitrary location which you want to find.

2.6368285629644522e-08 [-4.36512562  6.10654201  6.58197312] 9.98331343484 9.983313434838061

sensor1 [  0.   0.  10.] 6.38585280874e-08
[ 8.66025404  0.         -5.        ] 1.10316871818e-07
[-4.33012702  7.5        -5.        ] 7.93935118587e-08
[-4.33012702 -7.5        -5.        ] 1.07588488256e-07

The code that is used to pick a random time and position (this is not part of the solution)

```python
thetaPhi=np.array([[0,0],[2/3*np.pi,0],[2/3*np.pi,2/3*np.pi],[2/3*np.pi,4/3*np.pi]])
l=10 #m, distance from centre of apparatus to light sensor
p0=np.array([l*np.sin(thetaPhi[:,0])*np.cos(thetaPhi[:,1]),
    l*np.sin(thetaPhi[:,0])*np.sin(thetaPhi[:,1]),
    l*np.cos(thetaPhi[:,0])]).transpose()
print (p0)
print(np.sqrt(3)/2)

import random
t0=random.uniform(0,30e-9)

r0=10*random.uniform(0,1)**(1/3)
ct=random.uniform(-1,1)
phi=random.uniform(0,np.pi*2)
st=np.sqrt(1-ct*ct)

xpos=np.array([st*np.cos(phi),st*np.sin(phi),ct])*r0
print(t0,xpos,np.sqrt(np.dot(xpos,xpos)),r0)

vgroup=2.2e8

for p in p0:
    print(p, np.sqrt(np.dot(p-xpos,p-xpos))/vgroup+t0)
```

    [[  0.           0.          10.        ]
     [  8.66025404   0.          -5.        ]
     [ -4.33012702   7.5         -5.        ]
     [ -4.33012702  -7.5         -5.        ]]
    0.866025403784
    2.6368285629644522e-08 [-4.36512562  6.10654201  6.58197312] 9.98331343484 9.983313434838061
    [  0.   0.  10.] 6.38585280874e-08
    [ 8.66025404  0.         -5.        ] 1.10316871818e-07
    [-4.33012702  7.5        -5.        ] 7.93935118587e-08
    [-4.33012702 -7.5        -5.        ] 1.07588488256e-07

The equations we need to solve are
$$\frac{|\vec{s_i}-\vec{x_o}|}{v_g}=t-t_o$$

$$s_i^2+x_o^2+2\vec{s_i}\cdot\vec{x_o}=v_g^2\left(t^2-2tt_o+t_o^2\right)$$

# Solution for question 2

```python
import numpy as np
s=np.array(([[  0.,0.,10.,6.3859E-08],[8.66025404, 0., -5.,1.1032E-07],
             [-4.33012702, 7.5, -5.,7.9394E-08],[-4.33012702,-7.5, -5.,1.0759E-07]]))

vgroup=2.2e8

#using minimizer
def fun(par):
    xo=par[0:3]
    to=par[3]
    d2=0
    for sii in s:
        si=sii[0:3]
        ti=sii[3]
        d2+=(np.sqrt(np.dot(si-xo,si-xo))-vgroup*(ti-to))**2
    return d2

import scipy.optimize
x0=np.array([0,0,0,0])
sol=scipy.optimize.minimize(fun,x0,method='Powell')

# or using solver
def fun2(par):
    xo=par[0:3]
    to=par[3]
    d2=[]
    for sii in s:
        si=sii[0:3]
        ti=sii[3]
        d2.append(np.sqrt(np.dot(si-xo,si-xo))-vgroup*(ti-to))
    return np.array(d2)

print(s.shape)
sol2=scipy.optimize.fsolve(fun2,x0)
print('Solution from minimizer:',sol.x)
print('Solution from fsolve:',sol2, 'Check the solutions:',fun2(sol2))
```

    (4, 4)
    Solution from minimizer: [ -4.36595483e+00   6.10704000e+00   6.58259003e+00   2.63662486e-08]
    Solution from fsolve: [ -4.36595483e+00   6.10704000e+00   6.58259003e+00   2.63662486e-08] Check the solutions: [  3.76410014e-12  -1.25918831e-10  -1.64845915e-12   1.37077905e-10]

## Trajectories

An alpha-particle is a helium nucleus, with mass 6.644×10−27 kg or 4.002 u, and charge equal to twice the charge of an electron (but positive).  Calculate the trajectory of a 5MeV=(5*1.609 e-13 J) alpha particle as it moves by a gold nucleus (mass 196.966 u), with charge 79e, as a function of  the impact parameter b, below.  Assume both the alpha and the gold are point particles, and ignore special relativity.  Plot the scattering angle $\theta$ and energy loss of the alpha as a function of $b$, for values of b between 1e-16 and 1e-9m.

The force between two charge particles is given by the Coulomb potential:
$$\vec{F}=\frac{1}{4\pi\epsilon_0}\frac{q_1q_2 (\vec{r_2}-\vec{r_1})}{|\vec{r_2}-\vec{r_1}|^3}$$ with $\epsilon_0=8.85\times10^{-12}\frac{\rm{C}^2}{\rm{N\cdot m^2}}$

```python
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
q1=2*1.6e-19 # Coulombs, alpha
q2=79*1.6e-19 # coulombs, gold
k=1/(4*np.pi*8.85e-12)*q1*q2
mHe=6.644e-27
mAu=196.966/4.002*mHe
factor=(mAu+mHe)/mAu  # to transform x to gold-He distance
def traj(x,t):
    x0=factor*x[0:3]
    v0=x[3:6]
    distance=np.sqrt(np.dot(x0,x0))
    acceleration=k*x0/distance**3/mHe
    return np.array((v0,acceleration)).reshape(6)

E0=5*1.609e-13 #J
velocity=np.sqrt(2*E0/mHe)
init=np.array([-1.e-8,1e-14,0,velocity,0,0])
t0=np.linspace(0,-2*init[0]/velocity,500)
sol0=odeint(traj,init,t0)

blist=np.logspace(-16,-9,100)
sa=[]
for b in blist:
    init[1]=b
    sol=odeint(traj,init,t0,rtol=1e-6)
    angle=np.arctan2(sol[-1,4],sol[-1,3])*180/np.pi
    sa.append(angle)

plt.figure()
plt.title('Scattering angle versus impact parameter')
plt.ylabel('Scattering angle [degrees]')
plt.xlabel('Impact parameter(m)')
plt.semilogx(blist,sa)
plt.show()

```

# Don't forget to upload your work to eclass!

```python

```


text:
```
import numpy as np
```

# Lecture 2: Elimination with Matrices
講義動画：https://www.youtube.com/watch?v=QVKj3LADCnA

## 消去法 (Elimination)

$$
\begin{align}
x+2y+z &=2 \tag{1} \\
3x+8y+z&=12 \tag{2} \\
4y+z&=2 \tag{3}
\end{align}
$$

これを$Ax=b$の形にすると、$A$は次のように表される。
$$
A = \left[
  \begin{array}{ccc}
    1 & 2 & 1 \\
    3 & 8 & 1 \\
    0 & 4 & 1 \\
  \end{array}
\right]
$$

まず$(2)-(1)\times{}3$で$A_{2,1}$を0にする。→これで1列目のピボットができる。
$$
\left[
  \begin{array}{ccc}
    \bf{1} & 2 & 1 \\
    0 & 2 & -2 \\
    0 & 4 & 1 \\
  \end{array}
\right]
$$

次に、$(3)-(2)\times{}2$で$A_{3,2}$を0にする。→これで2列目と3列目のピボットができる。→これが上三角形(upper triangle)
$$
\left[
  \begin{array}{ccc}
  \bf{1} & 2 & 1 \\
  0 & \bf{2} & -2 \\
  0 & 0 & \bf{5} \\
  \end{array}
\right] = U
$$

ピボットの位置がすでに0だったら、それより下の行と入れ替える。

上記に$b$の列を追加してみる。

$$
\left[
  \begin{array}{ccc|c}
    1 & 2 & 1 & 2 \\
    3 & 8 & 1 & 12 \\
    0 & 4 & 1 & 2 \\
  \end{array}
\right]
$$

$(2)-(1)\times{}3$
$$
\left[
  \begin{array}{ccc|c}
    1 & 2 & 1 & 2 \\
    0 & 2 & -2 & 6 \\
    0 & 4 & 1 & 2 \\
  \end{array}
\right]
$$

$(3)-(2)\times{}2$
$$
\left[
  \begin{array}{ccc|c}
  1 & 2 & 1 & 2 \\
  0 & 2 & -2 & 6 \\
  0 & 0 & 5 & -10 \\
  \end{array}
\right]
$$

これにより、次の式を得る。
$$
\begin{align}
x+2y+z &=2 \\
2y-2z&=6 \\
5z&=-10
\end{align}
$$

ここから、$z=-2, y=1, x=2$が導き出せる。

行列×列＝列

$$
\left[
  \begin{array}{ccc}
    A_{1, 1} & A_{1, 2} & A_{1, 3} \\
    A_{2, 1} & A_{2, 2} & A_{2, 3} \\
    A_{3, 1} & A_{3, 2} & A_{3, 3} \\
  \end{array}
\right] \left[
  \begin{array}{c} a \\ b \\ c \end{array}
\right] = a\left[
  \begin{array}{c} A_{1,1} \\ A_{2,1} \\ A_{3,1} \end{array}
\right]+b\left[
  \begin{array}{c} A_{1,2} \\ A_{2,2} \\ A_{3,2} \end{array}
\right]+c\left[
  \begin{array}{c} A_{1,3} \\ A_{2,3} \\ A_{3,3} \end{array}
\right]
$$

行×行列＝行

$$
 \left[
  \begin{array}{ccc}
    a & b & c
  \end{array}
\right] \left[
  \begin{array}{ccc}
    A_{1, 1} & A_{1, 2} & A_{1, 3} \\
    A_{2, 1} & A_{2, 2} & A_{2, 3} \\
    A_{3, 1} & A_{3, 2} & A_{3, 3} \\
  \end{array}
\right] = a\left[
  \begin{array}{ccc} A_{1,1} & A_{2,1} & A_{3,1} \end{array}
\right]+b\left[
  \begin{array}{ccc} A_{2,1} & A_{2,2} & A_{2,3} \end{array}
\right]+c\left[
  \begin{array}{ccc} A_{3,1} & A_{3,2} & A_{3,3} \end{array}
\right]
$$

## 行列を使った消去法

上記の最初の変形（$\left[
  \begin{array}{ccc}
    1 & 2 & 1 \\
    3 & 8 & 1 \\
    0 & 4 & 1 \\
  \end{array}
\right]$
を
$\left[
  \begin{array}{ccc}
    1 & 2 & 1 \\
    0 & 2 & -2 \\
    0 & 4 & 1 \\
  \end{array}
\right]$
に変換）を行うには、行列の左側から
$\left[
  \begin{array}{ccc}
    1 & 0 & 0 \\
    -3 & 1 & 0 \\
    0 & 0 & 1 \\
  \end{array}
\right]$
を掛ける。

$$
\left[
  \begin{array}{ccc}
    1 & 0 & 0 \\
    -3 & 1 & 0 \\
    0 & 0 & 1 \\
  \end{array}
\right] \left[
  \begin{array}{ccc}
    1 & 2 & 1 \\
    3 & 8 & 1 \\
    0 & 4 & 1 \\
  \end{array}
\right] = \left[
  \begin{array}{ccc}
    1 & 2 & 1 \\
    0 & 2 & -2 \\
    0 & 4 & 1 \\
  \end{array}
\right]
$$

```
E21 = np.array([[1, 0, 0], [-3, 1, 0], [0, 0, 1]])
A = np.array([[1, 2, 1], [3, 8, 1], [0, 4, 1]])
A_ = np.dot(E21, A)
A_
```

    array([[ 1,  2,  1],
           [ 0,  2, -2],
           [ 0,  4,  1]])

さらに次の変形（$\left[
  \begin{array}{ccc}
    1 & 2 & 1 \\
    0 & 2 & -2 \\
    0 & 4 & 1 \\
  \end{array}
\right]$
を
$\left[
  \begin{array}{ccc}
    1 & 2 & 1 \\
    0 & 2 & -2 \\
    0 & 0 & 5 \\
  \end{array}
\right]$
に変換）を行うには、行列の左側から
$\left[
  \begin{array}{ccc}
    1 & 0 & 0 \\
    0 & 1 & 0 \\
    0 & -2 & 1 \\
  \end{array}
\right]$
を掛ける。

$$
\left[
  \begin{array}{ccc}
    1 & 0 & 0 \\
    0 & 1 & 0 \\
    0 & -2 & 1 \\
  \end{array}
\right] \left[
  \begin{array}{ccc}
    1 & 2 & 1 \\
    0 & 2 & -2 \\
    0 & 4 & 1 \\
  \end{array}
\right] = \left[
  \begin{array}{ccc}
    1 & 2 & 1 \\
    0 & 2 & -2 \\
    0 & 0 & 5 \\
  \end{array}
\right]
$$

```
E32 = np.array([[1, 0, 0], [0, 1, 0], [0, -2, 1]])
np.dot(E32, A_)
```

    array([[ 1,  2,  1],
           [ 0,  2, -2],
           [ 0,  0,  5]])

まとめると、$E_{32}(E_{21}A)=(E_{32}E_{21})A=U$

```
np.dot(np.dot(E32, E21), A)
```

    array([[ 1,  2,  1],
           [ 0,  2, -2],
           [ 0,  0,  5]])

## 置換行列 (Permutation Matrix)

左から行列を掛けると行の操作になる。
$$
\left[
  \begin{array}{cc}
    0 & 1 \\
    1 & 0 \\
  \end{array}
\right] \left[
  \begin{array}{cc}
    a & b \\
    c & d \\
  \end{array}
\right] = \left[
  \begin{array}{cc}
    c & d \\
    a & b \\
  \end{array}
\right]
$$

右から行列を掛けると列の操作になる。
$$
\left[
  \begin{array}{cc}
    a & b \\
    c & d \\
  \end{array}
\right] \left[
  \begin{array}{cc}
    0 & 1 \\
    1 & 0 \\
  \end{array}
\right] = \left[
  \begin{array}{cc}
    b & a \\
    d & c \\
  \end{array}
\right]
$$

## 逆行列
行列の変形をもとに戻す(undo)行列

$E=\left[
  \begin{array}{ccc}
    1 & 0 & 0 \\
    -3 & 1 & 0 \\
    0 & 0 & 1 \\
  \end{array}
\right]$
（2行目から1行目の3倍を引く）をもとに戻す行列は
$E^{-1}=\left[
  \begin{array}{ccc}
    1 & 0 & 0 \\
    3 & 1 & 0 \\
    0 & 0 & 1 \\
  \end{array}
\right]$

行列とその逆行列を掛けると単位行列になる。
$EE^{-1}=I$

```
E = np.array([[1, 0, 0], [-3, 1, 0], [0, 0, 1]])
E_inv = np.linalg.inv(E)
np.dot(E, E_inv)
```

    array([[1., 0., 0.],
           [0., 1., 0.],
           [0., 0., 1.]])


text:
# Mapping water extent and rainfall using Sentinel-1 and CHIRPS

* **Products used:**
[s1_rtc](https://explorer.digitalearth.africa/products/s1_rtc),
[rainfall_chirps_monthly](https://explorer.digitalearth.africa/products/rainfall_chirps_monthly)
**Keywords**: :index:`data used; sentinel-1`, :index:`data used; CHIRPS`, :index:`water; extent`, :index:`analysis; time series`
## Background

The United Nations have prescribed 17 "Sustainable Development Goals" (SDGs). This notebook attempts to monitor SDG Indicator 6.6.1 - change in the extent of water-related ecosystems. Indicator 6.6.1 has 4 sub-indicators:

    i. The spatial extent of water-related ecosystems
    ii. The quantity of water contained within these ecosystems
    iii. The quality of water within these ecosystems
    iv. The health or state of these ecosystems

This notebook primarily focuses on the first sub-indicator - spatial extents.

## Description

The notebook loads Sentinel-1 data, and using the sentinel-water-index, maps the spatial extent of water bodies. It also loads and plots monthly total rainfall from CHIRPS.  The last section will compare the water extent between two periods to allow visulazing where change is occuring.

***

## Load packages
Import Python packages that are used for the analysis.

```python
%matplotlib inline

import datacube
import xarray as xr
import matplotlib.pyplot as plt

from deafrica_tools.datahandling import load_ard
from deafrica_tools.dask import create_local_dask_cluster
from deafrica_tools.datahandling import wofs_fuser

from long_term_water_extent import (
    load_vector_file,
    get_resampled_labels,
    resample_water_observations,
    resample_rainfall_observations,
    calculate_change_in_extent,
    compare_extent_and_rainfall,
)
```

## Set up a Dask cluster

Dask can be used to better manage memory use and conduct the analysis in parallel.

```python
create_local_dask_cluster()
```

Client

Scheduler: tcp://127.0.0.1:38437

Dashboard: /user/chad/proxy/8787/status

Cluster

Workers: 1

Cores: 4

Memory: 28.14 GB

## Connect to Data Cube

```python
dc = datacube.Datacube(app="long_term_water_extent")
```

## Analysis parameters

The following cell sets the parameters, which define the area of interest and the length of time to conduct the analysis over.

* Upload a vector file for your water extent and your catchment to the `data` folder.
* Set the time range you want to use.
* Set the resampling strategy. Possible options include:
    *  `"1Y"` - Annual resampling, use this option for longer term monitoring
    * `"QS-DEC"` - Quarterly resampling from December
    * `"3M"` - Three-monthly resampling
    * `"1M"` - Monthly resampling

For more details on resampling timeframes, see the [xarray](https://xarray.pydata.org/en/v0.8.2/generated/xarray.Dataset.resample.html#r29) and [pandas](https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html#dateoffset-objects) documentation.

```python
water_extent_vector_file = "data/lake_baringo_extent.geojson"

water_catchment_vector_file = "data/lake_baringo_catchment.geojson"

time_range = ("2018-07", "2021") #earliest date of S1 is July 1st 2018

resample_strategy = "Q-DEC"

dask_chunks = dict(x=1000, y=1000)
```

## Get waterbody and catchment geometries

The next cell will extract the waterbody and catchment geometries from the supplied vector files, which will be used to load Water Observations from Space and the CHIRPS rainfall products.

```python
extent, extent_geometry = load_vector_file(water_extent_vector_file)
catchment, catchment_geometry = load_vector_file(water_catchment_vector_file)
```

## Load Sentinel-1 for Waterbody

The first step is to load the Sentinel-1 data using the extent geometry.

```python
extent_query = {
    "time": time_range,
    "resolution": (-20, 20),
    "output_crs": "EPSG:6933",
    "geopolygon": extent_geometry,
    "group_by": "solar_day",
    "dask_chunks":dask_chunks
}

ds = load_ard(dc=dc, products=["s1_rtc"], measurements=["vv", "vh"], **extent_query)
```

    Using pixel quality parameters for Sentinel 1
    Finding datasets
        s1_rtc
    Applying pixel quality/cloud mask
    Returning 196 time steps as a dask array

## Convert the Sentinel-1 digital numbers to dB

While Sentinel-1 backscatter is provided as linear intensiy, it is often useful to convert the backscatter to decible (dB) for analysis.
Backscatter in dB unit has a more symmetric noise profile and less skewed value distribution for easier statistical evaluation.

The Sentinel-1 backscatter data is converted from digital number (DN) to backscatter in decibel unit (dB) using the function:

\begin{equation}
10 * \log_{10}(\text{DN})
\end{equation}

```python
# Convert DN to db values.
ds["vv_db"] = 10 * xr.ufuncs.log10(ds.vv)
ds["vh_db"] = 10 * xr.ufuncs.log10(ds.vh)
```

## Calculate the SWI water index

The Sentinel-1A water index (SWI) is calculated as follows:

\begin{equation}
\text{SWI} = 0.1747 * \beta _{vv} + 0.0082 * \beta _{vh} * \beta _{vv} + 0.0023 * \beta _{vv}^{2} - 0.0015 * \beta _{vh}^{2} + 0.1904
\end{equation}

where  βvh and βvv represent the backscattering coefficient in VH polarization and VV polarization, respectively ([Tian et al., 2017](https://doi.org/10.3390/rs9060521)).

```python
# Calculate the Sentinel-1A Water Index (SWI).
ds['swi'] = (
    (0.1747 * ds.vv_db)
    + (0.0082 * ds.vh_db * ds.vv_db)
    + (0.0023 * ds.vv_db ** 2)
    - (0.0015 * ds.vh_db ** 2)
    + 0.1904
)

swi = ds[['swi']]
print(swi)
```

    <xarray.Dataset>
    Dimensions:      (time: 196, y: 1818, x: 852)
    Coordinates:
      * time         (time) datetime64[ns] 2018-07-01T15:56:48.569171 ... 2021-12...
      * y            (y) float64 9.433e+04 9.431e+04 ... 5.801e+04 5.799e+04
      * x            (x) float64 3.473e+06 3.473e+06 3.473e+06 ... 3.49e+06 3.49e+06
        spatial_ref  int32 6933
    Data variables:
        swi          (time, y, x) float32 dask.array<chunksize=(1, 1000, 852), meta=np.ndarray>
    Attributes:
        crs:           EPSG:6933
        grid_mapping:  spatial_ref

### Identify water in each resampling period

The second step is to resample the observations to get a consistent measure of the waterbody, and then calculate the classified as water for each period.

```python
resampled_water_ds, resampled_water_area_ds = resample_water_observations(
    swi, resample_strategy, radar=True
)

date_range_labels = get_resampled_labels(ds, resample_strategy)
```

    CPLReleaseMutex: Error = 1 (Operation not permitted)

    Automatic SWI threshold: 0.59

### Plot the change in water area over time

```python
fig, ax = plt.subplots(figsize=(15, 5))

ax.plot(
    date_range_labels,
    resampled_water_area_ds.values,
    color="red",
    marker="^",
    markersize=4,
    linewidth=1,
)
plt.xticks(date_range_labels, rotation=65)
plt.title(f"Observed Area of Water from {time_range[0]} to {time_range[1]}")
plt.ylabel("Waterbody area (km$^2$)")
plt.tight_layout()
```

## Load CHIRPS monthly rainfall

```python
catchment_query = {
    "time": time_range,
    "resolution": (-5000, 5000),
    "output_crs": "EPSG:6933",
    "geopolygon": catchment_geometry,
    "group_by": "solar_day",
    "dask_chunks":dask_chunks
}

rainfall_ds = dc.load(product="rainfall_chirps_monthly", **catchment_query)
```

### Resample to estimate rainfall for each time period

This is done by taking calculating the average rainfall over the extent of the catchment, then summing these averages over the resampling period to estimate the total rainfall for the catchment.

```python
catchment_rainfall_resampled_ds = resample_rainfall_observations(
    rainfall_ds, resample_strategy, catchment
)
```

## Compare waterbody area to catchment rainfall

This step plots the summed average rainfall for the catchment area over each period as a histogram, overlaid with the waterbody area calculated previously.

```python
figure = compare_extent_and_rainfall(
    resampled_water_area_ds, catchment_rainfall_resampled_ds, "mm", date_range_labels
)
```

### Save the figure

```python
figure.savefig("waterarea_and_rainfall_radar.png", bbox_inches="tight")
```

## Compare water extent for two different periods

For the next step, enter a baseline date, and an analysis date to construct a plot showing where water appeared, as well as disappeared, by comparing the two dates.

```python
baseline_time = "2018-07-01"
analysis_time = "2021-10-01"
```

```python
figure = calculate_change_in_extent(baseline_time, analysis_time, resampled_water_ds, radar=True)
```

### Save figure

```python
figure.savefig("waterarea_change_radar.png", bbox_inches="tight")
```

---

## Additional information

**License:** The code in this notebook is licensed under the [Apache License, Version 2.0](https://www.apache.org/licenses/LICENSE-2.0).
Digital Earth Africa data is licensed under the [Creative Commons by Attribution 4.0](https://creativecommons.org/licenses/by/4.0/) license.

**Contact:** If you need assistance, please post a question on the [Open Data Cube Slack channel](http://slack.opendatacube.org/) or on the [GIS Stack Exchange](https://gis.stackexchange.com/questions/ask?tags=open-data-cube) using the `open-data-cube` tag (you can view previously asked questions [here](https://gis.stackexchange.com/questions/tagged/open-data-cube)).
If you would like to report an issue with this notebook, you can file one on [Github](https://github.com/digitalearthafrica/deafrica-sandbox-notebooks).

**Compatible datacube version:**

```python
print(datacube.__version__)
```

    1.8.6

**Last Tested:**

```python
from datetime import datetime
datetime.today().strftime('%Y-%m-%d')
```

    '2022-01-10'

max_stars_repo_name: vikineema/deafrica-sandbox-notebooks

max_issues_repo_name: vikineema/deafrica-sandbox-notebooks

max_forks_repo_name: vikineema/deafrica-sandbox-notebooks