text:
# Clustering Techniques Writeup
### September 26, 2016

### K-Means Clustering

K-Means clustering is a regression technique which involves 'fitting' a number $n$ of given values from a dataset around a pre-defined number of $k$ clusters.  The K-means clustering process seeks to minimize the Euclidean distance from each point $n$ from its associated $k$ cluster.  Note that although the $k$ clusters inhabit the same coordinate space as the $n$ data points, the $k$ cluster locations are not necessarily derived from the set of $n$ points.

The process of K-means involves first defining the number of $k$ clusters.  Defining the number of $k$ means is a non-trivial problem that is dependent on the specific experiment at hand.  Analysis with an elbow plot showing explained variance as a function of number of clusers allows the designer to designate the $k$ number by choosing the point where adding additional $k$ clusters has diminishing returns for explained variance.

The actual algorithm behind K-means clustering is relatively simple.  Given $n$ data points and $k$ clusters, we're looking for µk (cluster means that minimize the Euclidean distance between the $n$ data points and said cluster).  Effectively, the objective function we're trying to minimize can be written as (courtesy SciKit):

$$
\begin{align}
\sum\limits_{i=0}^n \min_{µ_j \subseteq C} ||(x_j - µ_i)||^2
\end{align}
$$

For the analysis, we must first choose $k$ starting mean µ's.  This is often achieved by randomly sampling $k$ points from the dataset  arbitrarily.  Below, we have an image of the data set at the start, with arbitrarily placed means. *(Image sequence courtesy David Runyan, from Project Rhea)*

After choosing the initial set of starting means, we can apply a three step iterative process to eventually 'converge' on finalized means:

1)  **Fit n points to µ means:**   Calculate the Euclidean distances and assign each $n$ to its closest µ.  In this case, all points are closest to the blue node, so all $n$ are assigned to the blue µ.

2)  **Define new µ's by finding average:**   Find the center of each cluster and assign those as new µ's.  Only the blue µ moves, as the other nodes have no values assigned to them.

3)  **Repeat until convergence:**   Repeat until there is no change in the µ positions.  I've shown the final converged image below.

**PROS:**
1.  Oft-implemented, existing documentation.
2.  Has SciKit implementation in Python.

**CONS:**
1.  Defining number of means is a non-trivial process.
2.  Does not necessarily converge to global solution, as randomized starting positions will alter final result for different results.
3.  Does not set the means to positions to some value from the dataset (if we're interested in brain clusters, we want our clusters to be centered around some set of pre-defined nodes, not some weighted average determined by k-means.)
4.  Assumes similarly sized 'clusters', which is a very large assumption to be making for brain areas.
5.  Fails for anisotropically sorted data.

**Pseudocode for K-Means:**
**Algorithm 1:** $k_{means}$ - designate and find $k$ means in a dataset
**Input:** $k$ number of clusters, $n$ in x, y, z Euclidean-space data, $T$ tolerance for step size.
**Output:** x, y, z coordinates of $k_{means}$
1.
**procedure** $k_{means}(k, n, T)$
2.  $\ \ \ \ \mu_i =$ random_select$(k, n) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $ Randomly select k elements from n and assign them to set µ.
3.  $\ \ \ \ \mu_i = \min_{µ_j \subseteq C} ||(x_j - µ_i)||^2 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  \ \ \ \ \ \ \ \ \ \ \ \ \ \ $  Assign each $n$ to $\mu_i$ to the minimum values.
4.  $\ \ \ \ \mu_i = \sum\limits_{i=0}^n \min_{µ_j \subseteq C} average(n_i \ in \ C)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $ Move $\mu$ to average of x, y, z for all n in c.
5.  $\ \ \ $ when $\Delta$ in $\mu_i$ between steps 3/4 is $T$, return all $\mu_i \ \ \ \ \ \ $Repeat until $T$ tolerance change between average and assign.

### K-Medoid Clustering

K-Medoids clustering is similar to K-means in that both start with a known, pre-defined number of 'k' clusters.  However, K-medoids sets the 'k' clusters to points that occur naturally in the dataset of 'n' points.  Thus, while K-means finds new 'means' and fits each 'n' point to these new means by reducing some summed error, K-medoids instead seeks to maximize some measure of similarity (eg: the similarity matrix) between each points and its medoid (which itself is one of the pre-existing points).

The process for K-Medoids is similar to K-Means; again, the problem of defining the number of 'k' means is a difficult process.  From there, an additional difficulty lies in how to define the 'similarity' between two points.  A way to do so would be to make a similarity matrix that is has each corresponding value as the inverse of the Euclidean distance between the points.  Aside from this, a similar process would be applied iteratively to converge to medoids that maximize the similarity between their respective nodes.

The specific name of the algorithmic approach is Partitioning Around Medoids, or PAM.  The steps are, specifically:
1)  Arbitrarily select 'k' of the 'n' nodes as the medoids.
2)  Calculate the total 'similarity' (eg, by finding the inverse of the total of the distances between all 'n' and their closest 'k', or by using some other measure).
3)  Swap one 'n' with one 'k', and recalculate the 'similarity' measure.  If the 'similarity' increased, keep the new configuration and continue.  Otherwise, return to the previous configuration.

**PROS:**
1.  Simple, like K-means.

**CONS:**
1.  Has similiar faults as the K-means methodology (defining 'k', not necessarily global, assumes equal sized regions)
2.  Recalculating similarity at each step relative to all other points is very computationally intensive

**Pseudocode for K-Medoids:**
**Algorithm 2:** $k_{medoids}$ - designate and find $k$ medoids in a dataset
**Input:** $k$ number of clusters, $n$ in x, y, z Euclidean-space data
**Output:** x, y, z coordinates of $k_{medoids}$
1.  **procedure** $k_{medoids}(k, n, T)$
2.  $\ \ \ \ \mu_i =$ random_select$(k, n) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $ Randomly select k elements from n and assign them to set µ.
3.  $\ \ \ $ Similarity = distance($x_j - µ_i$) $\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  \ \ \ \ \ \ \ \ \ \ $  Calculate 'similarity' value by finding distance between point and $\mu$
4.  $\ \ \ $ Swap($\mu_i, \ x_j$) $\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $ Set $x_j$ to $\mu_i$.
5.  $\ \ \ $ when $\Delta$ Similarity is $T$, repeat for $\mu_{i + 1}$ and $x_j + 1 \ \ \ \ \ \ \ \ $Repeat until $T$ tolerance change between average and assign.

### Spectral Clustering
Up until now, we've focused on clustering techniques centered around compactness.  Spectral clustering is a different dataset fitting technique that focuses on connectivity instead.

The actual algorithm behind spectral clustering is somewhat more complicated. First and foremost, spectral clustering seeks to define some affinity matrix $A$.  Effectively, this affinity matrix grants a measure of the 'similarity' between two points.  Often times, this is accomplished by measuring the Euclidean distance between two points $x_i - x_j$ and using some inverse exponential (such as $exp(-\alpha \ ||x_i - x_j||^2)$, or $exp(-\alpha \ ||x_i - x_j||^2)/(2\sigma^2)$) to find $A_{i, \ j}$  Often times, $A_{i, \ j}$ is further limited by saying that for values where $||x_i - x_j||^2) < R$, then $A_{i, \ j} = 0$.

Now that we have this affinity matrix $A$, we want to find the Laplacian for the data. There are multiple different types of Laplacians; each is used after a specific normalization technique to find the diagonals to the $A_{i, \ j}$ matrix.  For this example, I'll focus on the simple Laplacian, which is just the degree matrix $D$ (a measure of the degree at a specific node) less the $A$ value.

Now that we have this Laplacian, the Laplacian can be broken down into blocks.  These subsets each have some individual eigenvalue/eigenvector, which can be used find these new clusters by linking all the points to their corresponding eigenvalue/eigenvectors.

**Pros:**
1. Spectral clustering techniques are capable of finding clusters that are near each other, but not distributed about some weighted center.
2. Does not necessarily need the (large) assumption that each cluster is  with some Gaussian distribution
2. Has a SciKit Python implementation.

**Cons:**
1. Still a bit confused as to how to calculate the intermediaries for the steps.

Sources:
Ng's paper (2012): http://ai.stanford.edu/~ang/papers/nips01-spectral.pdf,
Charles Martin blog: https://charlesmartin14.wordpress.com/2012/10/09/spectral-clustering/

**Pseudocode for Spectral Clustering:**
**Algorithm 3:** $spectral$ - designate and find $k$ means in a dataset
**Input:** $k$ number of clusters, $n$ in x, y, z Euclidean-space data
**Output:** x, y, z coordinates of $spectral$
1.  **procedure** $spectral(k, n, T)$
2.  $\ \ \ \ \mu_i =$ random_select$(k, n) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $ Randomly select k elements from n and assign them to set µ.
3.  $\ \ \ $ Similarity = distance($x_j, \ µ_i$) $\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  \ \ \ \ \ \ \ \ \ \ $  Calculate 'similarity' value as a function of distance, and assign as adjacency matrix
4.  $\ \ \ $ D = degrees(A) $\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $ Find degrees of similarity (adjacency) matrix A.
5.  $\ \ \ $ Find sparse matrix representation of D. $\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $ Ranks are the eigenvectors/values of the clusters.
6.  $\ \ \ $ Perform $k_{means}$ on eigenvalues/vectors $\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $Repeat until $T$ tolerance change between average and assign.

### Louvain Method of Community Detection
The Louvain method of community detection is a clustering technique that relies on monitoring the relative connectiveness of communities.

From Blondel, Guillaume, Lambiotte and Lefebvre's seminal paper that defined the Louvain method, they defined first a modularity value.  From the paper, "the modularity of a partition is a scalar value between -1 and 1 that measures the density of links inside communities as compared to links between communities."  In the case of weighted nodes (where for instance they used the number of calls between between two phone users), the modularity is defined as:

$$
\begin{align}
Q = \frac{1}{2m} \sum\limits_{i, \ j} \big[ A_{i, \ j} - \frac{k_i k_j}{2m} \big] \ \delta \big( c_i , c_j \big)
\end{align}
$$

In this case, the $A_{i, \ j}$ refers to the weight of the edge between $i$ and $j$.  From there, they defined $k_i$ as the $\sum_{j} A_{i, \ j}$ as the "weights of the edges attached to vertex $i$".  $c_i$ is the 'community' to which node $i$ is assigned, while the $\delta$-function $\delta \big( u, v \big)$ is 1 if  $u$ is $v$, and 0 otherwise.  $m$ is defined to be $\frac{1}{2}\sum_{i, \ j} A_{i, \ j}$  In lay speech, then, what this calculates is 'how connected' some node $i$ is relative some node $j$ (since we're subtracting from the edge weight of $ij$ the edge weights of all connections from $i$ and $j$).

From there, the paper progresses to explain how this traditional definition of modularity takes a significant amount of storage power to process.  Data analytics on the original data using continuous modularity calls take more processing power than is truly necessary.  What Blondel, Guillaume, Lambiotte, and Lefebvre did was define an algorithmically (and computationally) more efficient method to measure the *change* in modularity.  It does so by calculating a $\Delta Q$ value that is the *change* in modularity by moving (inserting) some node $i$ from an 'isolated neighborhood' into $C$:

$$
\begin{align}
\Delta Q  = \bigg[ \frac{\sum_{in} + \ k_{i, \ in}}{2m} - \big( \frac{\sum_{tot} + \ k_i}{2m} \big)^2 \bigg] - \bigg[ \frac{\sum_{in}}{2m} - \big( \frac{\sum_{tot}}{2m} \big)^2 - \big( \frac{k_i}{2m} \big)^2 \bigg]
\end{align}
$$

A similar equation is provided for the removal of a node $i$ from the neighborhood $C$.

Application-wise, the process for Louvain's method is similar to the approach for all the other regression techniques (K-means, K-medoids).  The approach first assigns and puts each node $N$ into a unique community (so it's in its own community).  They then calculate the $\Delta Q$ by adding the node $i$ into community $j$; if after all similar additions into other neighborhoods, the $\Delta Q$ is maximized by adding i to j, that process is fulfilled.  Otherwise, the system is reverted to before the addition.

By applying this until a local minimum occurs (when no changes to $i$ increase or decrease $\Delta Q$, the Louvain method can create a very solid mapping.

**Pros:**
1.  Computationally efficient, their calculations in the paper show a significant reduction in computation time (running on 6.3M nodes, their algorithm only took 197 seconds.  The only comparable algorithm that even converged was the Wakita and Tsumuri's CNN implementation method, which took almost three times as long).
2.  Accuracy seems to be high, according to their own discussion and results.

**Cons:**
1.  Requires weighted edges (we haven't generated these yet).  If our weights are defined by the number of connections within the epsilon ball, what's the point of using weighted edges over K-means or K-medoids?  Isn't that just an extra layer of complexity?
2.  Requires a C-based workaround that allows us to run the methods in Python (so I won't be able to use these lovely Markdown/Python notebooks).  Not too big a problem, but it's nice when packets like SciKit already has implementations of the other methods.  Link:  http://perso.crans.org/aynaud/communities/

Link to Louvain's Method paper: http://arxiv.org/abs/0803.0476

**Pseudocode for Louvain:**
**Algorithm 4:** $Louvain$ - designate and find $k$ means using Louvains in a dataset
**Input:** $k$ number of clusters, $n$ in x, y, z Euclidean-space data
**Output:** x, y, z coordinates of $Louvain$
1.  **procedure** $k_{medoids}(k, n, T)$
2.  $\ \ \ \ \mu_i =$ random_select$(k, n) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $ Randomly select k elements from n and assign them to set µ.
3.  $\ \ \ $ Similarity = distance($x_j - µ_i$) $\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  \ \ \ \ \ \ \ \ \ \ $  Calculate 'similarity' value by finding distance between point and $\mu$
4.  $\ \ \ $ Swap($\mu_i, \ x_j$) $\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $ Set $x_j$ to $\mu_i$ and measure change in $\Delta Q$.
5.  $\ \ \ $ when $\Delta Q$ decreases to $T$, repeat for $\mu_{i + 1}$ and $x_j + 1 \ \ \ \ \ \ $Repeat until $T$ tolerance in $\Delta Q$ change after adding/removing point from community.

### Table for Comparisons (summarizing each independent section)

    Clustering Type
    Pro
    Con

    K-means
    Oft-implemented, existing documentation.  Has SciKit implementation in Python.
    Defining number of means is a non-trivial process.  Does not necessarily converge to global solution, as randomized starting positions will alter final result for different results.  Does not set the means to positions to some value from the dataset (if we're interested in brain clusters, we want our clusters to be centered around some set of pre-defined nodes, not some weighted average determined by k-means.)  Assumes similarly sized 'clusters', which is a very large assumption to be making for brain areas.  Fails for anisotropically sorted data.

    K medoids
    Spectral clustering techniques are capable of finding clusters that are near each other, but not distributed about some weighted center.  Does not necessarily need the (large) assumption that each cluster is with some Gaussian distribution.  Has a SciKit Python implementation.
    Has similiar faults as the K-means methodology (defining 'k', not necessarily global, assumes equal sized regions).  Recalculating similarity at each step relative to all other points is very computationally intensive

    Spectral
    Spectral clustering techniques are capable of finding clusters that are near each other, but not distributed about some weighted center.  Does not necessarily need the (large) assumption that each cluster is with some Gaussian distribution  Has a SciKit Python implementation.
    Still a bit confused as to how to calculate the intermediaries for the steps.

    Louvain
    Computationally efficient, their calculations in the paper show a significant reduction in computation time (running on 6.3M nodes, their algorithm only took 197 seconds. The only comparable algorithm that even converged was the Wakita and Tsumuri's CNN implementation method, which took almost three times as long).
Accuracy seems to be high, according to their own discussion and results.
    Requires weighted edges (we haven't generated these yet). If our weights are defined by the number of connections within the epsilon ball, what's the point of using weighted edges over K-means or K-medoids? Isn't that just an extra layer of complexity?
Requires a C-based workaround that allows us to run the methods in Python (so I won't be able to use these lovely Markdown/Python notebooks). Not too big a problem, but it's nice when packets like SciKit already has implementations of the other methods.

```python

```

max_forks_repo_path: Jupyter/Clustering Techniques.ipynb

max_forks_repo_name: NeuroDataDesign/seelviz-archive


text:
# Chapter 6 Calculus

```python
from sympy import *
x, y, z = symbols('x y z')
init_printing(use_unicode=True)
```

## 6.1 微分

微分するには

```python
diff(cos(x),x)
```

```python
diff(exp(x**2),x)
```

とする。これらの関数はその場で微分操作を評価するので、未知関数には使えない。(微分方程式については後述)。

2階以上の微分をするには

```python
diff(x**4, x, x, x)
```

もしくは

```python
diff(x**4, x, 3) #xの3階微分.
```

とすればよい.

予想されるように、微分不可能になると

```python
diff(x**4, x, 5)
```

と0を返す.

多変数関数、たとえば$$e^{xyz}$$に対しても同様に

```python
expr = exp(x*y*z)
```

```python
diff(expr,x, y, y, z, z, z, z)
```

```python
diff(expr, x, y, 2, z, 4) #x,y:2階, z: 4階
```

```python
diff(expr, x, y, y, z, 4) #x, y, z: 4階
```

**コメント** `diff`は`Sympy`オブジェクトに対するメソッドとしても呼べる:

```python
expr.diff(x, y, y, z, 4)
```

#### 微分操作を行いたくないとき

上で述べたように、`diff`関数はその場で式を微分してしまうので、微分方程式には使えない.
その代わりに評価する前の微分を作るには`Derivative`を使えばよい。

```python
deriv = Derivative(expr, x, y, y, z, 4)
```

```python
deriv
```

これを評価するには`doit`メソッドを使う:

```python
deriv.doit() #diff関数で評価するものと同じ.
```

**コメント**: `Derivative()`は微分方程式などに使える.

## 6.2 積分

```python
integrate(cos(x), x)
```

**注意**: 不定積分であっても、未定定数は出力されない.

**例**: 定積分$$\int_{0}^{\infty}e^{-x} \, dx$$は

```python
integrate(exp(-x), (x, 0, oo)) #無限大は'o'２つ
```

と計算できる.

**例:** 2重積分$$\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}e^{-x^2-y^2}\, dx dy$$は

```python
integrate(exp(-x**2-y**2), (x, -oo, oo), (y, -oo, oo))
```

と計算できる.

#### 積分が評価できないとき

もし積分が評価できないときは, `Integral`で返す:

```python
expr = integrate(x**x, x)
```

```python
print(expr)
```

    Integral(x**x, x)

```python
expr
```

**コメント** 積分方程式に使える.

#### Integralの表式から積分値を評価したいとき

```python
expr = Integral(log(x)**2, x)
```

```python
expr
```

```python
expr.doit() #doit()をつければ評価される.
```

**コメント**: `integrate()`は合成関数とか、特殊関数に対しても使える!

```python
integ = Integral(x**y*exp(-x), (x, 0, oo))
```

```python
integ
```

```python
integ.doit()
```

収束性の場合分けまでやってくれる!!!

## 6.3 極限

極限$$\lim_{x \rightarrow 0}\frac{\sin(x)}{x}$$は

```python
limit(sin(x)/x, x, 0)
```

と評価できる.

```python
expr = x**2 / exp(x)
```

```python
expr.subs(x, oo) #oo/ooになる.
```

```python
limit(expr, x, oo)
```

*極限操作と代入操作は厳密に区別すべき！*

#### 極限を評価したくないとき

```python
expr = Limit((cos(x) - 1)/x, x, 0)
```

```python
expr
```

```python
expr.doit() #極限値を評価
```

#### 極限の方向を決めたいとき(正方向, 負方向)

$$\lim_{x\rightarrow 0^{+}}\frac{1}{x}$$

```python
limit(1/x, x, 0, '+') #第四引数に'+'のオプション
```

$$\lim_{x \rightarrow 0^{-}}\frac{1}{x}$$

```python
limit(1/x, x, 0, '-') #第四引数に'-'のオプション
```

## 6.4 べき級数

関数`f(x)`を`x0`周りで`n`次まで展開: `f(x).series(x, x0, n)` (`x0`, `n`は省ける)

```python
expr = exp(sin(x))
```

```python
expr.series(x, 0, 4) #最高次はオーダー記号として出力
```

```python
x + x**3 + x**6 + O(x**4) #大文字のO(オー)でオーダー記号を表せる.
```

```python
x*O(1)
```

**コメント** オーダータームを除去するには

```python
expr.series(x, 0, 4).removeO() #こちらも大文字のO(オー)
```

とすればよい. なお、オーダータームは原点以外の点の周りの展開でも現れる.

```python
exp(x - 6).series(x, x0=6) #パラメータ引数?
```

## 6.5 有限差分

省略
http://docs.sympy.org/latest/tutorial/calculus.html#derivatives参照

次は[Chapter7 Solvers](https://hiroyuki827.github.io/SymPy_tutorial/Chapter7_Solvers.html)へ!

ext: ipynb


text:
# COMP 135 day05 Lab: Univariate Gaussian distributions

## Outline

* **Part 1: Univariate Gaussian Distribution Basics**
* * Learn to sample from a univariate Gaussian

* **Part 2: ML Estimators and biased vs unbiased**
* * Simulate finite sample from true model, then estimate from that sample. Repeat many times.

* **Part 3: Properties of the univariate Gaussian**
* * Visualize linear transformations of Gaussian samples
* * Visualize sums of Gaussian samples
* * Visualize products of Gaussian samples
* * Learn to evaluate PDF of a univariate Gaussian and plot it alongside a histogram

## Skills

* To draw samples:
* * [`scipy.stats.norm(loc=..., scale=...).rvs(size=(nrows,ncols), random_state=SEED)`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html)
* * OR [`np.random.normal(loc=..., scale=...)`](https://numpy.org/doc/stable/reference/random/generated/numpy.random.normal.html)
* To evaluate the PDF:
* * Evaluate pdf with `scipy.stats.norm(loc=..., scale=...).pdf(values)`
* * Evaluate logpdf with `logpdf` instead of `pdf`

## Takeaways

* Linear transforms of Gaussian variables are also Gaussian distributed
* The sum of Gaussian variables is also Gaussian distributed
* The product of Gaussian variables is *not* Gaussian distributed

```python
import numpy as np
import pandas as pd
import scipy.stats
```

```python
np.set_printoptions(precision=3, suppress=False)
```

```python
pd.options.display.float_format = '{:,.3g}'.format  # show 3 digits of precision
```

```python
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style("whitegrid")
sns.set_context("notebook", font_scale=1.25)
```

```python

```

# Part 1: Univariate Gaussian basics

## Univariate Gaussian background

#### Random Variable: Scalar data $x \in \mathbb{R}$

This is a scalar real value.

#### Parameters: Location $\mu$ and scale $\sigma$

* Location or "mean" is just any real value: $\mu \in \mathbb{R}$
* Scale or standard-deviation $\sigma$ must be a positive value: $\sigma > 0$

We can also call *square* of the  "scale" parameter the *variance*: $\sigma^2$

#### Formal definition of PDF:

\begin{align}
\text{NormPDF}\left( x | \mu, \sigma^2 \right)
    = c(\mu, \sigma) \cdot e^{-\frac{1}{2} \frac{1}{\sigma^2} (x - \mu)^2 }
\end{align}

with normalizing "constant" (a term constant wrt our random variable $x$):

$$
c(\mu, \sigma) = \frac{1}{(2\pi)^{1/2}} \frac{1}{\sigma}
$$

##  Sampling from the "standard" univariate Gaussian

It is common to think about the Gaussian with zero mean and variance one.
Call this the "standard" normal or "standard" Gaussian.

We can use `np.random.randn` to sample from the standard normal

Set our random state so things are reproducible.

```python
prng = np.random.RandomState(11)
```

Draw several samples from the *standard* Normal (*standard* means zero mean and unit variance)

```python
prng.randn()
```

```python
prng.randn()
```

```python
prng.randn()
```

Draw 1000 samples and plot the histogram

```python
N = 1000
x_N = prng.randn(N)
sns.histplot(x_N);
plt.xlim([-3.5, 3.5]);
```

Do it again.

```python
N = 1000
x_N = prng.randn(N)
sns.histplot(x_N);
plt.xlim([-3.5, 3.5]);
```

##  Sampling from the *general* univariate Gaussian

Here, we'll assume that the mean and variance of the Gaussian are given

```python
mu = 11.22
sigma = 3.0
```

```python
prng.normal(loc=mu, scale=sigma)
```

```python
prng.normal(loc=mu, scale=sigma)
```

## Exercise 1a: Write code to plot a histogram of samples with $\mu = -1, \sigma = 3$

```python
N = 10000
# TODO draw samples (hint: use prng.normal with provided loc and scale)
# TODO plot histogram (hint: see lines above for calls to sns.histplot)
```

## Exercise 1b: Write code to plot a histogram of samples with $\mu = -5, \sigma = 0.2$

```python
N = 10000
# TODO draw samples (hint: use prng.normal with provided loc and scale)
# TODO plot histogram (hint: see lines above for calls to sns.histplot)
```

## Discussion 1c: Can you use the plots above to come up with reasonable rules for the "high density" regions of the univariate normal?

* A "super-majority" of samples (say 2/3 or ~67%) will occur between $\mu - a$ and $\mu + a$
* Almost all samples (say ~99%) will occur between $\mu - b$ and $\mu + b$

How should we pick $a$ and $b$? Can they be constants? Or should they depend on $\mu$ or $\sigma$?

## Exercise 1d: Try out your rules

```python
mu = 101.0
sigma = 15.0
a = 0.0 # TODO fix me
b = 0.0 # TODO fix me

N = 10000
x_N = prng.normal(loc=mu, scale=sigma, size=N)

n_match_for_a = np.sum(np.logical_and(x_N >= mu - a, x_N <= mu + a))
n_match_for_b = np.sum(np.logical_and(x_N >= mu - b, x_N <= mu + b))

print("%.2f samples occur between mu - %.2f and mu + %.2f" % (n_match_for_a/N, a, a))
print("%.2f samples occur between mu - %.2f and mu + %.2f" % (n_match_for_b/N, b, b))
```

```python

```

```python

```

# Part 2: Unbiased vs. Biased estimation of parameters

In pre-recorded videos, we developed the following estimators using the *maximum likelihood* principle:

\begin{align}
\hat{\mu} &= \frac{1}{N} \sum_{n=1}^N x_n
\\
\hat{\sigma}^2 &= \frac{1}{N} \sum_{n=1}^N (x_n - \hat{\mu} )^2
\end{align}

We further learned that:

* the ML-estimator for the mean is *unbiased*
* the ML-estimator for the variance is *biased*

Here, we're gonna work on *demonstrating* this bias (or lack there of) via simulation.

## Exercise 2a: Implement the ML estimation formulas in code

```python
def calc_ml_estimates_for_mean_and_variance(x_N):
    ''' Compute the ML estimators for Gaussian parameters given N observed data

    Args
    ----
    x_N : 1D array, shape (n_examples,) = (N,)
        N iid observations to be modeled as draws from a univariate Gaussian

    Returns
    -------
    mu_hat : float
    sigma_squared_hat : float, must be positive
    '''
    mu_hat = 0.0 # TODO fix me
    sigma_squared_hat = 1.0 # TODO fix me
    return mu_hat, sigma_squared_hat
```

## Exercise 2b: Implement simulations that repeatedly generate data from true model and then estimate parameters

```python
def run_many_simulations(n_trials=1000, N=10, mu_true=1.0, sigma_squared_true=1.0):
    ''' Produce estimates of parameters for many trials of a toy data experiment

    Each trial will:
    * Draw N samples of x_n ~ Normal(mean=mu_true, variance=sigma_squared_true)
    * Compute ML estimates using these samples

    Returns
    -------
    mu_hat_R : 1D array, shape (n_trials,) = (R,)
        Contains mu_hat estimate for each trial
    sigma_squared_hat_R : 1D array, shape (n_trials,) = (R,)
        Contains sigma_squared_hat estimate for each trial
    '''
    mu_hat_list = list()
    sigma_squared_hat_list = list()
    for trial in range(n_trials):
        x_N = np.zeros(N) # TODO draw N samples of x from a Gaussian with mu_true and sigma_squared_true
        assert x_N.shape == (N,)

        # Use your function from 2a to get estimates of the parameters
        mu_hat, sigma_squared_hat = calc_ml_estimates_for_mean_and_variance(x_N)

        # Record these estimates in a big list
        mu_hat_list.append(mu_hat)
        sigma_squared_hat_list.append(sigma_squared_hat)
    return np.asarray(mu_hat_list), np.asarray(sigma_squared_hat_list)

```

## Exercise 2c: Run simulations with R=200 trials and N=10 samples

We'll visualize the *distribution* of our estimator across trials.

We'll also look at how well the *average result* of these trials matches our *true* parameter.

Run the code below. No need to edit it.

```python
mu_true = 30.0;             # Keep this at 30
sigma_squared_true = 10.0;  # Keep this at 10
N = 10
R = 200

mu_hat_R, sigma_squared_hat_R = run_many_simulations(
    n_trials=R, N=N, mu_true=mu_true, sigma_squared_true=sigma_squared_true)
```

### Visualize true $\mu$ and estimates from many trials

Just run the code below to produce a visual. You don't need to edit it.

```python
mu_hat_mean = np.mean(mu_hat_R);
print("Empirical average value of \hat{\mu} = %.2f over %d trials" % (mu_hat_mean, R))

# Plot histogram of mu_hats across all R trials
sns.histplot(mu_hat_R, color='g', alpha=0.4);

# Plot "true" mu and "mean" mu as vertical lines
plt.plot(mu_hat_mean * np.ones(2), np.asarray([0, plt.gca().get_ylim()[1]]), 'g--', label='empirical $E[\hat{\mu}]$');
plt.plot(mu_true * np.ones(2), np.asarray([0, plt.gca().get_ylim()[1]]), 'b--', label='$\mu_{true}$');
plt.legend(bbox_to_anchor=(1.0, 0.8));

plt.title("Distribution of $\hat{\mu}$ over R=%d trials with N=%d" % (R, N));
```

### Compare true $\sigma^2$ and estimates from many trials

Just run the code below to produce a visual. You don't need to edit it.

```python
sigma_squared_hat_mean = np.mean(sigma_squared_hat_R);
print("Empirical average value of \hat{\sigma}^2 = %.2f" % sigma_squared_hat_mean)

sns.histplot(sigma_squared_hat_R, color='m', alpha=0.2);

# Plot "true" sigma and "mean" sigma_hat as vertical lines
plt.plot(sigma_squared_hat_mean * np.ones(2), np.asarray([0, plt.gca().get_ylim()[1]]), 'm--', label='$E[ \hat{\sigma}^2]$');
plt.plot(sigma_squared_true * np.ones(2), np.asarray([0, plt.gca().get_ylim()[1]]), 'r--', label='$\sigma^2_{true}$');
plt.legend(loc='upper right');

plt.title("Distribution of $\hat{\sigma}^2$ over R=%d trials with N=%d" % (R, N));
```

## Exercise 2d: Repeat the above with $R=10000$ trials (still keep $N=10$)

TODO copy code here and edit as needed.

## Discussion 2e: Do you see the bias for sigma? Do you see a "lack" of bias for mu?

TODO discuss. What visually in the plots above reveals "bias"? Do the patterns you see change with many trials?

## Discussion 2f: Does the amount of bias agree with derivations?

Remember, we suggested that:

$$
\mathbb{E}[ \hat{\sigma}^2 ] = \frac{N-1}{N} \sigma_{\text{true}}
$$

TODO discuss

```python

```

# Part 3: Visualizing transformations of Gaussians

### Exercise 3a: Linear transformations of Gaussians

* 1) Drawn 10000 samples of $X$ a standard Normal.
* 2) Transform each $X$ into $Y$, using linear transform $Y \gets 3 X + 2$
* 3) Draw the histograms of the resulting distributions
* 4a) Add a PDF plot for ML-estimated parameters using the samples of $X$
* 4b) Add a PDF plot for ML-estimated parameters using the samples of $Y$

Note that there are coding TODOs only for 1 and 4

```python
# STEP 1: Draw samples
x_N = prng.randn(10000)

# STEP 2: Apply linear transform: Y = 3 * X + 2
y_N = prng.rand(10000) # TODO replace with transform

# STEP 3: Draw histograms
fig, axgrid = plt.subplots(nrows=2, ncols=1, sharex=True, sharey=True, figsize=(10, 5))
sns.histplot(x_N, color='b', alpha=0.2, ax=axgrid[0], bins=np.linspace(-8, 12, 51), stat='density');
sns.histplot(y_N, color='r', alpha=0.2, ax=axgrid[1], bins=np.linspace(-8, 12, 51), stat='density');

# STEP 4a: Compute ML-estimates of parameters for X
mu_hat, sigma_squared_hat = calc_ml_estimates_for_mean_and_variance(x_N)
xgrid_G = np.linspace(-8, 12, 301);
pdfgrid_G = np.zeros(xgrid_G.size) # TODO call scipy.stats.norm's pdf function to get the pdf
axgrid[0].plot(xgrid_G, pdfgrid_G, 'b-', label="ML fit: $\mu$=%.2f $\sigma^2$=%.2f" % (mu_hat, sigma_squared_hat));
axgrid[0].legend(loc='upper right');
axgrid[0].set_title('$X \sim \mathcal{N}(0, 1)$')

# STEP 4b: Compute ML-estimates of parameters for Y
mu_hat, sigma_squared_hat = calc_ml_estimates_for_mean_and_variance(y_N)
ygrid_G = np.linspace(-8, 12, 301);
pdfgrid_G = np.zeros(ygrid_G.size) # TODO call scipy.stats.norm's pdf function to get the pdf
axgrid[1].plot(ygrid_G, pdfgrid_G, 'r-', label="ML fit: $\mu$=%.2f $\sigma^2$=%.2f" % (mu_hat, sigma_squared_hat));
axgrid[1].legend(loc='upper right');
axgrid[1].set_title('$Y = 3 X + 2$');

plt.tight_layout();
```

## Discussion 3b: Does Y appear to be Gaussian distributed?

TODO discuss. What clues did you use?

```python

```

## Exercise 3c: Sums of Gaussian random variables

* 1a) Draw 10000 samples of $X$ from a standard Normal.
* 1b) Draw 10000 samples of $Y$ from a standard Normal.

* 2) Apply a sum transform: $Z = X+ Y$
* 3) Draw the histograms of the resulting distributions

* 4) Add a PDF plot for ML-estimated parameters using the samples of $Z$

Note that there are coding TODOs only for 1 and 4

```python
# STEP 1: Draw samples
x_N = prng.randn(10000)
y_N = prng.randn(10000)

# STEP 2: Apply sum transform: Z = X + Y
z_N = prng.rand(10000) # TODO replace with transform

# STEP 3: Draw histograms
xmin = -4
xmax = +4
fig, axgrid = plt.subplots(nrows=3, ncols=1, sharex=True, sharey=True, figsize=(10, 8))
sns.histplot(x_N, color='b', alpha=0.2, ax=axgrid[0], bins=np.linspace(xmin, xmax, 51), stat='density');
sns.histplot(y_N, color='y', alpha=0.2, ax=axgrid[1], bins=np.linspace(xmin, xmax, 51), stat='density');
sns.histplot(z_N, color='g', alpha=0.2, ax=axgrid[2], bins=np.linspace(xmin, xmax, 51), stat='density');

# STEP 4a: Compute ML-estimates of parameters for X
mu_hat, sigma_squared_hat = calc_ml_estimates_for_mean_and_variance(x_N)
xgrid_G = np.linspace(xmin, xmax, 301);
pdfgrid_G = np.zeros(ygrid_G.size) # TODO call scipy.stats.norm's pdf function to get the pdf
axgrid[0].plot(xgrid_G, pdfgrid_G, 'b-', label="ML fit: $\mu$=%.2f $\sigma^2$=%.2f" % (mu_hat, sigma_squared_hat));
axgrid[0].legend(loc='upper right');
axgrid[0].set_title('$X \sim \mathcal{N}(0, 1)$')

# STEP 4b: Compute ML-estimates of parameters for Y
mu_hat, sigma_squared_hat = calc_ml_estimates_for_mean_and_variance(y_N)
ygrid_G = np.linspace(xmin, xmax, 301);
pdfgrid_G = np.zeros(ygrid_G.size) # TODO call scipy.stats.norm's pdf function to get the pdf
axgrid[1].plot(ygrid_G, pdfgrid_G, 'y-', label="ML fit: $\mu$=%.2f $\sigma^2$=%.2f" % (mu_hat, sigma_squared_hat));
axgrid[1].legend(loc='upper right');
axgrid[1].set_title('$Y \sim \mathcal{N}(0, 1)$');

# STEP 4b: Compute ML-estimates of parameters for Z
mu_hat, sigma_squared_hat = calc_ml_estimates_for_mean_and_variance(z_N)
zgrid_G = np.linspace(xmin, xmax, 301);
pdfgrid_G = np.zeros(ygrid_G.size) # TODO call scipy.stats.norm's pdf function to get the pdf
axgrid[2].plot(zgrid_G, pdfgrid_G, 'g-', label="ML fit: $\mu$=%.2f $\sigma^2$=%.2f" % (mu_hat, sigma_squared_hat));
axgrid[2].legend(loc='upper right');
axgrid[2].set_title('$Z = X + Y$');
plt.tight_layout();
```

## Discussion 3d: Does Z = X + Y above appear to be Gaussian distributed?

```python
# TODO discuss!
```

## Exercise 3e: Products of Gaussian random variables

* 1a) Draw 10000 samples of $X$ from a standard Normal.
* 1b) Draw 10000 samples of $Y$ from a standard Normal.

* 2) Apply a *product* transform: $Z = X * Y$
* 3) Draw the histograms of the resulting distributions

* 4) Add a PDF plot for ML-estimated parameters using the samples of $Z$

```python
# STEP 1: Draw samples
x_N = prng.randn(10000)
y_N = prng.randn(10000)

# STEP 2: Apply product transform: Z = X * Y
z_N = prng.rand(10000) # TODO replace with transform

# STEP 3: Draw histograms
xmin = -3
xmax = +3
fig, axgrid = plt.subplots(nrows=3, ncols=1, sharex=True, sharey=True, figsize=(10, 8))
sns.histplot(x_N, color='b', alpha=0.2, ax=axgrid[0], bins=np.linspace(xmin, xmax, 51), stat='density');
sns.histplot(y_N, color='y', alpha=0.2, ax=axgrid[1], bins=np.linspace(xmin, xmax, 51), stat='density');
sns.histplot(z_N, color='g', alpha=0.2, ax=axgrid[2], bins=np.linspace(xmin, xmax, 51), stat='density');

# STEP 4a: Compute ML-estimates of parameters for X
mu_hat, sigma_squared_hat = calc_ml_estimates_for_mean_and_variance(x_N)
xgrid_G = np.linspace(xmin, xmax, 301);
pdfgrid_G = np.zeros(xgrid_G.size) # TODO compute PDF
axgrid[0].plot(xgrid_G, pdfgrid_G, 'b-', label="ML fit: $\mu$=%.2f $\sigma^2$=%.2f" % (mu_hat, sigma_squared_hat));
axgrid[0].legend(loc='upper right');
axgrid[0].set_title('$X \sim \mathcal{N}(0, 1)$')

# STEP 4b: Compute ML-estimates of parameters for Y
mu_hat, sigma_squared_hat = calc_ml_estimates_for_mean_and_variance(y_N)
ygrid_G = np.linspace(xmin, xmax, 301);
pdfgrid_G = np.zeros(xgrid_G.size) # TODO compute PDF
axgrid[1].plot(ygrid_G, pdfgrid_G, 'y-', label="ML fit: $\mu$=%.2f $\sigma^2$=%.2f" % (mu_hat, sigma_squared_hat));
axgrid[1].legend(loc='upper right');
axgrid[1].set_title('$Y \sim \mathcal{N}(0, 1)$');

# STEP 4b: Compute ML-estimates of parameters for Z
mu_hat, sigma_squared_hat = calc_ml_estimates_for_mean_and_variance(z_N)
zgrid_G = np.linspace(xmin, xmax, 301);
pdfgrid_G = np.zeros(xgrid_G.size) # TODO compute PDF
axgrid[2].plot(zgrid_G, pdfgrid_G, 'g-', label="ML fit: $\mu$=%.2f $\sigma^2$=%.2f" % (mu_hat, sigma_squared_hat));
axgrid[2].legend(loc='upper right');
axgrid[2].set_title('$Z = X * Y$');
plt.tight_layout();
```

## Discussion 3f: Does Z = X * Y above appear to be Gaussian distributed?

TODO discuss!

```python

```


text:
```python
import numpy as np
import pandas as pd
import linearsolve as ls
import matplotlib.pyplot as plt
plt.style.use('classic')
%matplotlib inline
```

# Class 18: Exogenous Shocks in the New-Keynesian Model

In this notebook, we will compute impulse responses to exogenous shocks in the new-Keynesian model and interpret the results. Then we compute a stochastic simulation of the model.

## The New-Keynesian Model

The most basic version of the New-Keynesian Model can be expressed as:

\begin{align}
y_t & = E_t y_{t+1}  - \left( r_{t} - \bar{r}\right) + g_t\\
i_{t} & = r_{t} + E_t \pi_{t+1}\\
i_{t} & = \bar{r} + \pi^T + \phi_{\pi}\big(\pi_t - \pi^T\big) + \phi_{y}\big(y_t - \bar{y}\big) + v_t\\
\pi_t -\pi^T & = \beta \left( E_t\pi_{t+1} - \pi^T\right)  + \kappa (y_t -\bar{y})+ u_t,
\end{align}

where: $y_t$ is (log) output, $r_t$ is the real interest rate, $i_t$ is the nominal interest rate, $\pi_t$ is the rate of inflation between periods $t-1$ and $t$, $\bar{r}$ is the long-run average real interest rate or the *natural rate of interest*, $\beta$ is the household's subjective discount factor, and $\pi^T$ is the central bank's inflation target. The coeffieints $\phi_{\pi}$ and $\phi_{y}$ reflect the degree of intensity to which the central bank *endogenously* adjusts the nominal interest rate in response to movements in inflation and output.

The variables $g_t$, $u_t$, and $v_t$ represent exogenous shocks to aggregate demand, inflation, and monetary policy. They follow AR(1) processes:

\begin{align}
g_{t+1} & = \rho_g g_{t} + \epsilon^g_{t+1}\\
u_{t+1} & = \rho_u u_{t} + \epsilon^u_{t+1}\\
v_{t+1} & = \rho_v v_{t} + \epsilon^v_{t+1}.
\end{align}

The goal is to compute impulse responses in the model to all three exogenous shocks. We will use the following parameterization:

| $$\bar{y}$$ | $$\beta$$ | $$\bar{r}$$    | $$\kappa$$ | $$\pi^T$$ | $$\phi_{\pi}$$ | $$\phi_y$$ | $$\rho_g$$ | $$\rho_u$$ | $$\rho_v$$ |
|-------------|-----------|----------------|------------|-----------|----------------|------------|------------|------------|------------|
| 0           | 0.995     | $$-\log\beta$$ | 0.25       | 0.02/4    | 1.5            | 0.5/4      | 0.5        | 0.5        | 0.5        |

## Solution

```python
# Create a variable called 'parameters' that stores the model parameter values in a Pandas Series. CELL PROVIDED
parameters = pd.Series(dtype=float)
parameters['y_bar'] = 0
parameters['beta'] = 0.995
parameters['r_bar'] = -np.log(parameters.beta)
parameters['kappa'] = 0.25
parameters['pi_T'] = 0.02/4
parameters['phi_pi'] = 1.5
parameters['phi_y'] = 0.5/4
parameters['rho_g'] = 0.5
parameters['rho_u'] = 0.5
parameters['rho_v'] = 0.5

# Create variable called 'var_names' that stores the variable names in a list with state variables ordered first
var_names = ['g','u','v','y','pi','i','r']

# Create variable called 'shock_names' that stores an exogenous shock name for each state variable.
shock_names = ['e_g','e_u','e_v']

# Define a function that evaluates the equilibrium conditions of the model solved for zero.
def equilibrium_equations(variables_forward,variables_current,parameters):

    # Parameters
    p = parameters

    # Current variables
    cur = variables_current

    # Forward variables
    fwd = variables_forward

    # IS equation
    is_equation = fwd.y - (cur.r -p.r_bar) + cur.g - cur.y

    # Fisher_equation
    fisher_equation = cur.r + fwd.pi - cur.i

    # Monetary policy
    monetary_policy = p.r_bar + p.pi_T + p.phi_pi*(cur.pi - p.pi_T) + p.phi_y*cur.y + cur.v - cur.i

    # Phillips curve
    phillips_curve = p.beta*(fwd.pi- p.pi_T) + p.kappa*cur.y + cur.u - (cur.pi-p.pi_T)

    # Demand process
    demand_process = p.rho_g*cur.g - fwd.g

    # Monetary policy process
    monetary_policy_process = p.rho_v*cur.v - fwd.v

    # Inflation process
    inflation_process = p.rho_u*cur.u - fwd.u

    # Stack equilibrium conditions into a numpy array
    return np.array([
            is_equation,
            fisher_equation,
            monetary_policy,
            phillips_curve,
            demand_process,
            monetary_policy_process,
            inflation_process
        ])

# Initialize the model into a variable named 'nk_model'
nk_model = ls.model(equations = equilibrium_equations,
                 n_states=3,
                 var_names=var_names,
                 shock_names=shock_names,
                 parameters = parameters)

# Compute the steady state numerically using .compute_ss() method of nk_model
guess = [0,0,0,0,0.01,0.01,0.01]
nk_model.compute_ss(guess)
```

```python
# Find the log-linear approximation around the non-stochastic steady state and solve using .approximate_and_solve() method of nk_model
# set argumement 'log_linear' to False
nk_model.approximate_and_solve(log_linear=False)

# Print the solution to the model
print(nk_model.solved())
```

    Solution to the linear system:

    g[t+1] = 0.5·g[t]+e_g[t+1]

    u[t+1] = 0.5·u[t]+e_u[t+1]

    v[t+1] = 0.5·v[t]+e_v[t+1]

      y[t] = 0.8909·g[t]-1.7729·u[t]-0.8909·v[t]

     pi[t] = 0.4432·g[t]+1.108·u[t]-0.4432·v[t]

      i[t] = 0.7762·g[t]+1.4404·u[t]+0.2238·v[t]

      r[t] = 0.5546·g[t]+0.8864·u[t]+0.4454·v[t]

**Questions**

1. Will a positive shock to $g_t$ increase or decrease $y_t$?
2. Will a positive shock to $g_t$ increase or decrease $\pi_t$?
3. Will a positive shock to $u_t$ increase or decrease $y_t$?
4. Will a positive shock to $u_t$ increase or decrease $\pi_t$?

**Answers**

1. The coefficient on $g_t$ in the solution for $y_t$ is *positive* so the shock will *increase* $y_t$
2. The coefficient on $g_t$ in the solution for $\pi_t$ is *positive* so the shock will *increase* $\pi_t$
1. The coefficient on $u_t$ in the solution for $y_t$ is *negative* so the shock will *decrease* $y_t$
4. The coefficient on $u_t$ in the solution for $\pi_t$ is *positive* so the shock will *increase* $\pi_t$

## Impulse Responses

Compute a 21 period impulse response of the model's variables to a:

* 0.01 unit shock to $g_t$
* 0.01/4 unit shock to $u_t$
* 0.01/4 unit shock to $v_t$

```python
# Compute impulse responses
nk_model.impulse(T=21,t0=5,shocks=[0.01,0.01/4,0.01/4])
```

### Demand Shock

Plot the computed impulses responses of the nominal interest rate, the real interest rate, output, and inflation to the 0.01 unit demand shock. Express inflation and interest rates in *annualized* (e.g., multiplied by 4) terms.

```python
# Create figure. CELL PROVIDED
fig = plt.figure(figsize=(12,8))

# Create upper-left axis
ax1 = fig.add_subplot(2,2,1)
# Create upper-right axis
ax2 = fig.add_subplot(2,2,2)
# Create lower-left axis
ax3 = fig.add_subplot(2,2,3)
# Create lower-right axis
ax4 = fig.add_subplot(2,2,4)

# Set axis 1 ylabel
ax1.set_ylabel('% dev from steady state')
# Set axis 2 ylabel
ax2.set_ylabel('% dev from steady state')
# Set axis 3 ylabel
ax3.set_ylabel('% dev from steady state')
# Set axis 4 ylabel
ax4.set_ylabel('% dev from steady state')

# Set axis 1 limits
ax1.set_ylim([-0.5,3.5])
# Set axis 2 limits
ax2.set_ylim([-0.5,3.5])
# Set axis 3 limits
ax3.set_ylim([-0.5,2])
# Set axis 4 limits
ax4.set_ylim([-0.5,2])

# Plot the nominal interest rate, real interest rate, output, and inflation
(nk_model.irs['e_g']['i']*400).plot(ax=ax1,lw=4,alpha=0.75,title='Nominal Interest (annualized)',grid=True)
(nk_model.irs['e_g']['r']*400).plot(ax=ax2,lw=4,alpha=0.75,title='Real Interest (annualized)',grid=True)
(nk_model.irs['e_g']['y']*100).plot(ax=ax3,lw=4,alpha=0.75,title='Output',grid=True)
(nk_model.irs['e_g']['pi']*400).plot(ax=ax4,lw=4,alpha=0.75,title='Inflation (annualized)',grid=True)
```

### Inflation Shock

Plot the computed impulses responses of the nominal interest rate, the real interest rate, output, and inflation to the 0.01/4 unit inflation shock. Express inflation and interest rates in *annualized* (e.g., multiplied by 4) terms.

```python
# Create figure. CELL PROVIDED
fig = plt.figure(figsize=(12,8))

# Create upper-left axis
ax1 = fig.add_subplot(2,2,1)
# Create upper-right axis
ax2 = fig.add_subplot(2,2,2)
# Create lower-left axis
ax3 = fig.add_subplot(2,2,3)
# Create lower-right axis
ax4 = fig.add_subplot(2,2,4)

# Set axis 1 ylabel
ax1.set_ylabel('% dev from steady state')
# Set axis 2 ylabel
ax2.set_ylabel('% dev from steady state')
# Set axis 3 ylabel
ax3.set_ylabel('% dev from steady state')
# Set axis 4 ylabel
ax4.set_ylabel('% dev from steady state')

# Set axis 1 limits
ax1.set_ylim([-.5,2])
# Set axis 2 limits
ax2.set_ylim([-0.5,2])
# Set axis 3 limits
ax3.set_ylim([-0.5,0.5])
# Set axis 4 limits
ax4.set_ylim([-0.5,1.5])

# Plot the nominal interest rate, real interest rate, output, and inflation
(nk_model.irs['e_u']['i']*400).plot(ax=ax1,lw=4,alpha=0.75,title='Nominal Interest (annualized)',grid=True)
(nk_model.irs['e_u']['r']*400).plot(ax=ax2,lw=4,alpha=0.75,title='Real Interest (annualized)',grid=True)
(nk_model.irs['e_u']['y']*100).plot(ax=ax3,lw=4,alpha=0.75,title='Output',grid=True)
(nk_model.irs['e_u']['pi']*400).plot(ax=ax4,lw=4,alpha=0.75,title='Inflation (annualized)',grid=True)
```

### Monetary Policy Shock

Plot the computed impulses responses of the nominal interest rate, the real interest rate, output, and inflation to the 0.01/4 unit monetary policy shock. Express inflation and interest rates in *annualized* (e.g., multiplied by 4) terms.

```python
# Create figure. CELL PROVIDED
fig = plt.figure(figsize=(12,8))

# Create upper-left axis
ax1 = fig.add_subplot(2,2,1)
# Create upper-right axis
ax2 = fig.add_subplot(2,2,2)
# Create lower-left axis
ax3 = fig.add_subplot(2,2,3)
# Create lower-right axis
ax4 = fig.add_subplot(2,2,4)

# Set axis 1 ylabel
ax1.set_ylabel('% dev from steady state')
# Set axis 2 ylabel
ax2.set_ylabel('% dev from steady state')
# Set axis 3 ylabel
ax3.set_ylabel('% dev from steady state')
# Set axis 4 ylabel
ax4.set_ylabel('% dev from steady state')

# Set axis 1 limits
ax1.set_ylim([-0.2,0.6])
# Set axis 2 limits
ax2.set_ylim([-0.2,0.6])
# Set axis 3 limits
ax3.set_ylim([-0.5,0.1])
# Set axis 4 limits
ax4.set_ylim([-0.5,0.1])

# Plot the nominal interest rate, real interest rate, output, and inflation
(nk_model.irs['e_v']['i']*400).plot(ax=ax1,lw=4,alpha=0.75,title='Nominal Interest (annualized)',grid=True)
(nk_model.irs['e_v']['r']*400).plot(ax=ax2,lw=4,alpha=0.75,title='Real Interest (annualized)',grid=True)
(nk_model.irs['e_v']['y']*100).plot(ax=ax3,lw=4,alpha=0.75,title='Output',grid=True)
(nk_model.irs['e_v']['pi']*400).plot(ax=ax4,lw=4,alpha=0.75,title='Inflation (annualized)',grid=True)
```

## Example: Stochastic Simulation

Compute a 401 period stochastic simulation of the model.

| $$\rho_g$$ | $$\sigma_g$$ | $$\rho_u$$ | $$\sigma_u$$ | $$\rho_v$$ | $$\sigma_v$$ |
|------------|--------------|------------|--------------|------------|--------------|
| 0.5        | 0.015        | 0.5        | 0.003        | 0.5        | 0.001        |

The parameter values are chosen so that the following simulated values are on par with the historical data for the US:

* standard deviation of output
* standard deviation of inflation
* standard deviation of the nominal interest rate
* correlation between output and inflation
* correlation between output and the nominal interest rate

```python
# Form the covariance matrix for the exogenous white noise processes
Sigma = [[0.015**2,0,0],[0,0.003**2,0],[0,0,0.001**2]]
```

```python
# Compute the stochastic simulation using the .stoch_sim() method of nk_model. Set seed to 126.
nk_model.stoch_sim(T=201,cov_mat=Sigma,seed=126)
```

```python
# Create 12x4 figure. PROVIDED
fig = plt.figure(figsize=(12,12))

# Create axis. PROVIDED
ax1 = fig.add_subplot(3,1,1)

# Plot the output gap
ax1.plot(nk_model.simulated['y']*100)
ax1.set_title('Output gap')
ax1.set_ylabel('Percent')
ax1.grid()

# Create axis. PROVIDED
ax2 = fig.add_subplot(3,1,2)

# Plot inflation (annualized)
ax2.plot(nk_model.simulated['pi']*400)
ax2.set_title('Inflation (annualized)')
ax2.set_ylabel('Percent')
ax2.grid()

# Create axis. PROVIDED
ax3 = fig.add_subplot(3,1,3)

# Plot the nominal interest rate (annualized)
ax3.plot(nk_model.simulated['i']*400)
ax3.set_title('Nominal interest rate (annualized)')
ax3.set_ylabel('Percent')
ax3.grid()
```

```python
# Compute standard deviations of output (times 100), inflation (times 400), and the nominal interest rate (times 400)
nk_model.simulated[['y','pi','i']].std()*[100,400,400]
```

    y     1.686375
    pi    3.779509
    i     6.242343
    dtype: float64

```python
# Compute correlation coefficients of output, inflation, and the nominal interest rate
nk_model.simulated[['y','pi','i']].corr()
```

      y
      pi
      i

      y
      1.000000
      0.697158
      0.756250

      pi
      0.697158
      1.000000
      0.993517

      i
      0.756250
      0.993517
      1.000000


text:
+ This notebook is part of lecture 18 *Properties of determinants* in the OCW MIT course 18.06 by Prof Gilbert Strang [1]
+ Created by me, Dr Juan H Klopper
    + Head of Acute Care Surgery
    + Groote Schuur Hospital
    + University Cape Town
    + Email me with your thoughts, comments, suggestions and corrections

Linear Algebra OCW MIT18.06 IPython notebook [2] study notes by Dr Juan H Klopper is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.

+ [1] OCW MIT 18.06
+ [2] Fernando Pérez, Brian E. Granger, IPython: A System for Interactive Scientific Computing, Computing in Science and Engineering, vol. 9, no. 3, pp. 21-29, May/June 2007, doi:10.1109/MCSE.2007.53. URL: http://ipython.org

```python
from IPython.core.display import HTML, Image
css_file = 'style.css'
HTML(open(css_file, 'r').read())
```

<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Philosopher:400,700,400italic,700italic' rel='stylesheet' type='text/css'>

```python
from sympy import init_printing, Matrix, symbols
from IPython.display import HTML
from warnings import filterwarnings
```

```python
init_printing(use_latex = 'mathjax')
filterwarnings('ignore')
```

```python
a, b, c = symbols('a b c')
```

# Properties of the determinant

## Notation

* The determinant of a matrix A is written as det(A) or |A|

## Main properties

* There are three main properties (first three listed below) and seven that follow from them

* det(I)=1
* A row exchange changes the sign of the determinant
    * We now know the determinant of every permutation matrix
* Multiplying any row with a constant results in the determinant also being multiplied by that constant
        * Only works when altering a single row (the determinant is a linear function of each row separately)
        * An alternate way of seeing this is
$$ \begin{vmatrix} a+a' & b+b' \\ c & d \end{vmatrix}=\begin{vmatrix} a & b \\ c & d \end{vmatrix}+\begin{vmatrix} a' & b' \\ c & d \end{vmatrix} $$
* If two rows are equal then the determinant is zero
    * This follows from property number two, where if we interchange rows the sign must change
    * This only works for zero, since row exchange leaves the matrix unchanged, which now can't have a different determinant (opposite sign)
* Subtracting a constant multiple of one row from another leaves the determinant unchanged
    * This flows from property three above
$$ \begin{vmatrix} a & b \\ c-la & d-lb \end{vmatrix}=\begin{vmatrix} a & b \\ c & d \end{vmatrix}+\begin{vmatrix} a & b \\ -la & -lb \end{vmatrix} $$
    * From property four (determinant of matrix with two similar rows equal zero) we now have the following
$$ =\begin{vmatrix} a & b \\ c & d \end{vmatrix}+\left( -l \right) \begin{vmatrix} a & b \\ a & b \end{vmatrix}\\ =\begin{vmatrix} a & b \\ c & d \end{vmatrix}+0\\ =\begin{vmatrix} a & b \\ c & d \end{vmatrix} $$
* The determinant of a matrix with a complete row (or columns) of zero is zero
    * This also follows from property three above, but multiplying a row by zero
* The determinant of an upper triangular matrix is the product of the elements of the main diagonal (the pivots)
    * Same goes for a diagonal matrix
        * This helps us to develop the expression for the determinant of a 2×2 matrix
$$ \begin{bmatrix} a & b \\ c & d \end{bmatrix}\\ \begin{bmatrix} \frac { c }{ a } a & \frac { c }{ a } b \\ c & d \end{bmatrix}\\ \begin{bmatrix} a & b \\ c-\frac { c }{ a } a & d-\frac { c }{ a } b \end{bmatrix}\\ \begin{bmatrix} a & b \\ 0 & d-\frac { c }{ a } b \end{bmatrix}\\ \therefore \quad \begin{vmatrix} a & b \\ 0 & d-\frac { c }{ a } b \end{vmatrix}\\ =\left( a \right) \left( d-\frac { c }{ a } b \right) \\ =ad-a\frac { c }{ a } b\\ =ad-bc $$
    * We can change an upper triangular matrix into a diagonal matrix by row operations (leaving the determinant unchanged by property five)
    * Now we can use the first part of the third property and take out each pivot
    * Eventually we are left with the identity matrix and the product of all the pivots
    * For a zero on the main diagonal we can use the property of a row of zeros and know the determinant is zero
* If the determinant is zero, the matrix is singular (only has the zero solution / not invertible)
* The determinant of the product of matrices
$$ \left| AB \right| =\left| A \right| \left| B \right| \\ \left| { A }^{ -1 } \right| =\frac { 1 }{ \left| A \right|  } \\ \left| { A }^{ 2 } \right| =\left| A \right| \left| A \right| ={ \left| A \right|  }^{ 2 }\\ \left| cA \right| ={ c }^{ n }\left| A \right|  $$
* For the determinant of a transpose of a matrix we have
$$ \left| { A }^{ T } \right| =\left| A \right|  $$

## Example problems

* Do the following by making use of the properties above

### Example problem 1

```python
A = Matrix([[101, 201, 301], [102, 202, 302], [103, 203, 303]])
A
```

$$\left[\begin{matrix}101 & 201 & 301\\102 & 202 & 302\\103 & 203 & 303\end{matrix}\right]$$

#### Solution

* By constant multiple subtraction we get
$$ \begin{bmatrix} 101 & 201 & 301 \\ 1 & 1 & 1 \\ 1 & 1 & 1 \end{bmatrix} $$
* Two identical rows, thus the determinant is zero

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

$$0$$

### Example problem 2

```python
A = Matrix([[1, a, a ** 2], [1, b, b ** 2], [1, c, c ** 2]])
A
```

$$\left[\begin{matrix}1 & a & a^{2}\\1 & b & b^{2}\\1 & c & c^{2}\end{matrix}\right]$$

#### Solution

* Subtracting constant multiple of row 1 from rows 2 and three and expanding the elements
$$ =\begin{vmatrix} 1 & a & { a }^{ 2 } \\ 0 & b-a & { b }^{ 2 }-{ a }^{ 2 } \\ 0 & c-a & { c }^{ 2 }-{ a }^{ 2 } \end{vmatrix}\\ =\begin{vmatrix} 1 & a & { a }^{ 2 } \\ 0 & b-a & \left( b-a \right) \left( b+a \right)  \\ 0 & c-a & \left( c-a \right) \left( c+a \right)  \end{vmatrix} $$
* Using property three that states that the determinant is a linear property of each row
$$ =\left( b-a \right) \left( c-a \right) =\begin{vmatrix} 1 & a & { a }^{ 2 } \\ 0 & 1 & \left( b+a \right)  \\ 0 & 1 & \left( c+a \right)  \end{vmatrix} $$
* Another elimination on row three
$$ \left( b-a \right) \left( c-a \right) =\begin{vmatrix} 1 & a & { a }^{ 2 } \\ 0 & 1 & \left( b+a \right)  \\ 0 & 1 & \left( c+a \right)  \end{vmatrix} $$
* Now we have upper triangular form and the determinant is the product of the elements in the main diagonal and also multiplying the (*b*-*a*)(*c*-*a*)
$$ =\left( b-a \right) \left( c-a \right) \left( b-c \right)  $$

```python
(A.det()).factor() # Calculating the determinant and factorizing the result
```

$$- \left(a - b\right) \left(a - c\right) \left(b - c\right)$$

* This is called a Vandermonde matrix
* http://en.wikipedia.org/wiki/Vandermonde_matrix

### Example problem 3

```python
A = Matrix([1, 2, 3]) * Matrix([[1, -4, 5]])
A
```

$$\left[\begin{matrix}1 & -4 & 5\\2 & -8 & 10\\3 & -12 & 15\end{matrix}\right]$$

#### Solution

* The rows of the resultant 3×3 matrix is linearly dependent, i.e. they are 1 times the row (1,-4,5), then twice this same row for row two and lastly three times the same row for row three
* This means that the determinant will be zero

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

$$0$$

### Example problem 4

```python
A = Matrix([[0, 1, 3], [-1, 0, 4], [-3, -4, 0]])
A
```

$$\left[\begin{matrix}0 & 1 & 3\\-1 & 0 & 4\\-3 & -4 & 0\end{matrix}\right]$$

#### Solution

* Note how this matrix is *skew* symmetric
* This means that AT=-A
* With the matrices AT and -A being equal, their determinant are equal
* Remember, though that the determinant of a matrix is the same as the determinant of the transpose of the matrix
$$ \left| A \right| =\left| { A }^{ T } \right| =\left| -A \right| ={ \left( -1 \right)  }^{ 3 }\left| A \right| =-\left| A \right| \\ \left| A \right| =-\left| A \right| \\ \therefore \quad \left| A \right| =0 $$

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

$$0$$

* Not all skew symmetric matrices have a zero determinant
* It only works because *n* is odd for this size matrix being 3×3 allowing for the negative

```python

```

max_stars_repo_name: solomonxie/jupyter-notebooks

max_issues_repo_name: solomonxie/jupyter-notebooks

max_forks_repo_name: solomonxie/jupyter-notebooks


text:
# Cálculo

## Gradiente

El gradiente de $f$ respecto a $A\in \mathbb{R}^{m\times n}$ es la matriz de derivadas parciales, de manera que

\begin{equation}
(\nabla_A f(A))_{ij}=\frac{\partial f(A)}{\partial A_{ij}}
\end{equation}

El tamaño de $\nabla_A f(A)$ es siempre del tamaño de $A$. El gradiente es una función únicamente definida si la función $f$ retorna un valor escalar. Además,

- $\nabla_{x}(f(x)+g(x))=\nabla_x f(x)+\nabla_x g(x)$

- Para $t\in\mathbb{R},\nabla_x (tf(x))=t\nabla_x f(x)$

Los gradientes se pueden ver como generalizaciones de las derivadas parciales a múltiples variables.

# La Hessiana.

Supongan que $f:\mathbb{R}^n \rightarrow \mathbb{R}$ es una función que toma un vector en $\mathbb{R}^n$ y regresa un número real. Entonces, la Hessiana es una matriz respecto a x, escrita $\nabla_x^2 f(x)$ o simplemente $H$ tal que

\begin{equation}
(\nabla_x^2 f(x))_{ij}=\frac{\partial^2 f(x)}{\partial x_i \partial x_j}
\end{equation}

Se puede notar que la Hessiana siempre es simétrica.

# Mínimos cuadrados

Suponiendo $A\in \mathbb{R}^{m\times n}$ y un vector $b\in R^m$ tal que $b\notin \mathcal{R}(A)$. En esta situación no es posible encontrar un vector $x\in\mathbb{R}^n$ tal que $Ax=b$, así que entonces, se quiere encontrar un vector $x$ tal que $Ax$ sea tan cercano como es posible a $b$, medido como la norma cuadrada euclidiana $||Ax-b||_2^2$.

Usando que $||x||_{2}^2=x^T x$, se tiene que

\begin{equation}
||Ax-b||_{2}^2=(Ax-b)^T(Ax-b)=x^TA^TAx-2b^TAx+b^Tb
\end{equation}

Tomando el gradiente con respecto a $x$ y usando las propiedades de los gradientes, se tiene que

\begin{equation}
\nabla_x (x^TA^TAx-2b^TAx+b^Tb)=2A^T Ax-2A^T b
\end{equation}

Haciendo esta expresión igual a cero y resolviendo para $x$, se tienen las ecuaciones normales.

\begin{equation}
x=(A^T A)^{-1}A^T b
\end{equation}

```python

```


text:
Machine Learning 101

## La necesidad de Encontrar Patrones
### Estimando la probabilidad del movimiento en el valor de activos

### Detectando _outliers_ en comportamiento de mercados

### Distribuciones Multimodales

### Creación de Nueva Información

la junta de gobierno esta a favor  de con otivo por las posi la pol itica mantenticar las a 2019 la inclativo del economia de conticiones a catarias aciver con extancialente se argendo de contisus inacticas delmedcano de resis nicresi pora contenticas inacimosteres respuntaria la potribal para contiyan la polisicam entrencialecienos

## ¿Qué es Machine Learning (*ML*)?
> *Machine Learning es el área de estudio que le da a las computadoras la habilidad de aprender sin necesidad de ser explicitamente programadas.* - Arthur Samuel

> _[Machine learning is] a set of methods that can automatically detect patterns in data, and then use the uncovered patterns to predict future data, or to perform other kinds of decision making under uncertainty._ -Kevin P. Murphy

```python
# Cargamos las librerías que estaremos ocupandos para la clase
import matplotlib.pyplot as plt
import numpy as np
from numpy.random import seed, randn
from ipywidgets import interact, IntSlider, FloatSlider
```

```python
%config InlineBackend.figure_format = "retina"
plt.rcParams["figure.figsize"] = (10, 5)
```

En general, se dice que un programa de computación aprende de una experiencia **E** con respecto a una _clase_ de tarea **T** y medida de desempeño **P** si el desempeño en **T**, medido bajo **P** mejora con experiencia **E**

## La tarea $T$

Machine learning nos ayuda a atacar problemas o tareas que serían difícil para un humano escribir usando métodos *tradicionales* de programación, i.e., definiendo la estructura lógica que debería seguir un programa para llegar al resultado.

Una **tarea** en machine learning es descrita en términos de como el sistema debe procesar una **observación** ($x$), la cuál definimos como un arreglo de números.

$$
    x = \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix},
$$

el cuál conocemos como un vector. Un vector, en python, es representado por medio de un numpy array.

```python
np.array([1, 2, 3, 5])
```

Si, por ejemplo, nuestras observaciones son las características de una persona, podríamos describir una observación de la siguiente manera:

$$
   \mbox{persona}_i = \begin{bmatrix} \mbox{edad} \\ \mbox{peso} \\ \mbox{estatura} \\ \mbox{color ojos} \\ \vdots \end{bmatrix}
$$

¿Cuáles podrían ser las tareas?

La **tarea**, dado $\mbox{persona}_i$, sería el resultado que pudiesemos obtener una vez procesado esta información: esperanza de vida de la persona; prima de un seguro; probabilidad de gusto hacia cierto producto, etc.

**Tareas en Machine Learning**
- Enseñar a un robot a tocar el piano
- Reconocer los dígitos de una imagen
- Detectar transacciones fradulentas
- Estimar ventas futuras

**Tareas no apropiadas para machine learning (aunque posibles)**
- Escribir un programa que calcule la sequencia fibonacci
- Escribir un programa que arroje reportes
- Escribir un programa de operaciones aritméticas

**Una tarea que no necesita machine learning.**

Supongamos que, dadas tres observaciones de personas (`persona1`, `persona2`, `persona3`), ¿sería la tarea **T** de cácular el índice de masa corporal un problema de Machine Learning?

$$
    \mbox{BMI} = \frac{\mbox{Masa}_{Kg}}{\mbox{Altura}_{m}^2}
$$

donde,
$$
    \mbox{persona}_i = \begin{bmatrix} \mbox{masa} \\ \mbox{altura} \end{bmatrix}
$$

```python
import numpy as np
persona1 = np.array([60, 1.40])
persona2 = np.array([180, 1.20])
persona3 = np.array([62, 1.66])
```

```python
# En otras palabras, la tarea para un algoritmo de ML es encontrar la función apropiada
# para los datos en cuestión. En este ejemplo no es necesario usar ML dado que
# estamos seguros del resultado que arroja la tarea dado cada ejemplo.

def bmi(persona):
    """
    Calcula el Indice de masa corporal de una persona

    Parametros
    -----------
    Persona: list
        lista con cuyo primer elemento es el peso (en Kg)
        y el segundo la altura (en metros) de una persona
    """
    return persona[0] / persona[1] ** 2

persona1 = np.array([60, 1.40])
persona2 = np.array([180, 1.20])
persona3 = np.array([62, 1.66])

print(bmi(persona1))
print(bmi(persona2))
print(bmi(persona3))
```

    30.612244897959187
    125.0
    22.49963710262738

### Las clases de Tareas:
- Clasificación
    - Asignar una calificación del 1 al 5 para una película
    - Sentimiento de mercado respecto a un producto
- Regresión
    - Estimar ventas futuras
    - Estimar precio de productos
    - Estimar peso de una persona (dada su altura)
- *Machine Translation*
    - Traducir
- Detección de anomalías
    - Dado un sistema, encontrar patrones estadísticos fuera del promedio de observaciones "normales".
- *Denoising*
    - Dada información corrompida, limpiar los datos
- Reducción de Dimensiones

## La Experiencia $E$

La *experiencia* es la información de la cual se alimenta nuestro modelo para realizar estimaciones. El tipo de experiencia depende de la clase de tarea de la cuál estemos hablando. Esta información se ve representada en un conjunto de observaciones que llamamos un **dataset** (la base de datos).

Considerando el ejemplo anterior, si contamos con 3 ejemplos de personas, el agrupamiento de cada una de las personas conformaría nuestra base de datos.

$$
    \mbox{personas} = \begin{bmatrix}
    \mbox{persona}_1, \mbox{persona}_2, \mbox{persona}_3
    \end{bmatrix}
$$

$$
    \mbox{personas} = \begin{bmatrix}
        \mbox{edad}_1 & \mbox{edad}_2 & \mbox{edad}_3 \\
        \mbox{peso}_1 & \mbox{peso}_2 & \mbox{peso}_3 \\
        \mbox{estatura}_1 & \mbox{estatura}_2 & \mbox{estatura}_3 \\
        \mbox{color ojos}_1 & \mbox{color ojos}_2 & \mbox{color ojos}_3 \\
        \vdots & \vdots & \vdots
    \end{bmatrix}
$$

Otra representación, más común, de los datos es representando cada columna las carácterísticas y, cada fila, las observaciones.

Persona | edad | peso | estatura
--- | --- | --- | ---
$\mbox{persona}_1$ | $\mbox{edad}_1$ | $\mbox{peso}_1$ | $\mbox{estatura}_1$
$\mbox{persona}_2$ | $\mbox{edad}_2$ | $\mbox{peso}_2$ | $\mbox{estatura}_2$
$\mbox{persona}_3$ | $\mbox{edad}_3$ | $\mbox{peso}_3$ | $\mbox{estatura}_3$

Dependiendo del tipo de experiencia que un modelo necesita categorizamos un modelo como **supervised** o **unsupervised** learning.

Hablando en términos generales,**unsupervised learning** se encarga de buscar patrones dentro de nuestro *dataset*. Para este tipo de problemas, no sabemos de antemano que es lo que agrupa diferentes grupos. En este tipo de problemas solo conocemos las características (o **features**) de nuestra base de datos. Por otro lado, un modelo de **supervised learning** contiene tanto características como un factor diferenciador: un precio, una clase, etc.

```python
from sklearn.datasets import make_moons

X, y = make_moons(n_samples=200, noise=0.1, random_state=1643)
fig = plt.figure(figsize=(15,5))
ax1, ax2 = fig.add_subplot(121), fig.add_subplot(122)
ax1.scatter(*X.T)
ax2.scatter(*X.T, c=y, cmap=plt.cm.viridis)
ax1.set_title("Unsupervised Learning")
ax2.set_title("Supervised Learning")
plt.show()
```

## El Desempeño $P$

Para evaluar las habilidades de un modelo en machine learning es necesaria una manera cuantitativa de medir su **desempeño**. Una medida de desempeño es tan importante como el problema a solucionar.

La finalidad de un desempeño se divide en dos:
1. $P$ nos debe ayudar a llegar a una $T$ óptima para la base de datos
2. $P$ mide como nuestro modelo, bajo $T$, generaliza las decisiones para ejemplos **aún no observados** (predicciones)

En general, dada una seríe de datos, dividimos nuestros datos en dos partes: **training** y **test**. El objetivo es entonces encontrar $T$ dado **training** que generalice los valores no observados dentro de test.

Fig. Partición de la base de datos entre training y test

### Medidas de desempeño:

#### Clasificación
$$
    \mbox{accuracy} = \frac{\#\mbox{Predicciones Correctas}}{\#\mbox{Valores totales}}
$$

#### Regresión

$$
    \mbox{MSE} = \frac{1}{n}\sum_{i=1}^n (\hat{y}_i - y_i)^2
$$

$$
    \mbox{R}^2 = 1 - \frac{\mbox{RSS}}{\mbox{TSS}} = 1 -  \frac{\sum_i^n(y_i - \hat y_i)^2}{\sum_i^n(y_i - \bar y)^2}
$$

Ejercicio

Supngamos que entrenado un modelo de clasificación binario ($\{0, 1\}$) nos arroja predicciones dentro de la lista `preds`. Si los valores verdaderos se encuentran dentro de la lista `verdaderos`, calcula el *accuracy* del modelo.

```python
predicciones = np.array([1, 1, 0, 1, 0, 0])
verdadero = np.array([1, 0, 0, 1, 1, 0])
```

## La Regresión Lineal (v.1)

Como su nombre menciona, la regresión lineal es un problema de regresión. En este modelo asumimos que nuestros **features** se relacionan con los **targets** de una manera lineal.

Otra manera de plantear el problema es decir que nuestras observaciones son las *variables independientes* y la variable a predecir la *variable dependiente*

Como un primer ejemplo, consideremos que tenemos una base de datos (**E**) $(x, y)$.

x|y
-|-
$x_1$|$y_1$
$x_2$|$y_2$
$x_3$|$y_3$
...|...
$x_{100}$|$y_{100}$

Nuestra hipótesis (la tarea $T$ a resolver), es poder aproximar $y$ dado $x$ de la siguiente manera:

$$
    \hat y = \theta_0 + \theta_1 x
$$

```python
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import seed, randn
# Semilla inicial para llegar al mismo resultado
seed(1643)
# Número de observaciones simuladas
nobs = 100
# Variables independientes
X = np.linspace(0, 1, nobs) + randn(nobs) * 0.3
# Variables dependientes
theta0, theta1 = 0.3, 1.5 # --- No los conocemos
y = theta0 + theta1 * X + randn(nobs) * 0.35
plt.scatter(X, y)
plt.xlabel("$X$")
plt.ylabel("$y$")
plt.show()
```

¿Cómo encontramos $\theta_0$, $\theta_1$?

**¡Necesitamos una medida de desempeño $P$!**

Consideramos $P$ como la función
$$
\begin{align}
    J(\theta_0, \theta_1) &= \frac{1}{n}\sum_{i=1}^{n} (\theta_0 + \theta_1 x_i - y_i)^2 \\
                          &= \frac{1}{n} \left[(\theta_0 + \theta_1 x_1 - y_1)^2 + (\theta_0 + \theta_1 x_2 - y_2)^2 + \ldots + (\theta_0 + \theta_1 x_{n} - y_{n})^2\right]
\end{align}
$$

Nos gustaría minimizar $J(\theta_0, \theta_1)$

Ejercicio

Si `theta0_est = 1.3` y `theta1_est = 0.3` son nuestras estimaciones de $\theta_0$ y $\theta_1$ dado `X` y `y`, cálcula $J($ `theta0_est`, `theta1_est`$)$

¿Por qué $J(\theta_0, \theta_1)$?

1. Mide las desviaciones al cuadrado entre nuestra hipótesis y la observación real
2. Minimiza el negativo de la log-verosimilitud

```python
def J(yhat, y): return np.sum((yhat - y) ** 2) / len(y)

def plotJ(theta0, theta1):
    """
    Función para graficar la norma L2 de una regresión lineal
    """
    seed(31415)
    x = np.array([1, 1.5, 2, 3, 4, 5, 5.6])
    xrange = np.linspace(min(x), max(x), 100)
    y = 1 + 1.5 * x + randn(len(x))
    yhat = theta0 + theta1 * xrange
    plt.scatter(x, y)
    plt.plot(xrange, yhat, color="crimson")
    plt.ylim(0, 12)
    plt.title(r"$J(\theta_0, \theta_1) =$" + format(J(theta0 + x * theta1 , y), "0.2f"))
    for xi, yi in zip(x, y):
        pred = theta0 + xi * theta1
        plt.vlines(xi, yi, pred)
        plt.text(xi + 0.1, yi, round((yi - pred) ** 2, 2))
```

```python
interact(plotJ,
        theta0=FloatSlider(min=-1, max=2, step=0.1, continuous_update=False),
        theta1=FloatSlider(min=0, max=3, step=0.1, continuous_update=False))
```

    interactive(children=(FloatSlider(value=0.0, continuous_update=False, description='theta0', max=2.0, min=-1.0)…

    <function __main__.plotJ(theta0, theta1)>

¿Por qué $J(\theta_0, \theta_1)$ tiene sentido?

##### Óptimo Global (Un único mínimo)

```python
t0, t1 = np.mgrid[-3:3:0.05, -5:5:0.05]
z = (t0.reshape(-1, 1) + t1.reshape(-1, 1) *  X - y) ** 2
z = z.sum(axis=1).reshape(t0.shape)
plt.contourf(t0, t1, z, 25)
plt.xlabel(r"$\theta_0$", fontsize=15)
plt.ylabel(r"$\theta_1$", fontsize=15)
plt.title(r"$J(\theta_0, \theta_1)$", fontsize=15)
plt.colorbar();
```

Para un *training dataset* $(x, y)$ con $n$ observaciones, se puede demostrar que los los valores $\theta_0^*$, $\theta_1^*$ que minimizan  $J(\theta_0, \theta_1)$ son:

$$
    \theta_1^* = \frac{n\sum_{i=1}^nx_iy_i - \left(\sum_{i=1}^nx_i\right)\left(\sum_{i=1}^ny_i\right)}{n\sum_{i=1}^nx_i^2 - \left(\sum_{i=1}^nx_i\right)^2}
$$

$$
    \theta_0^* = \frac{1}{n}\sum_{i=1}^ny_i - \theta_1^* \frac{1}{n}\sum_{i=1}^nx_i
$$

Ejercicio

Calcula `theta0_opt`, `theta1_opt` consierando las ecuaciones de arriba y compara los resultados contra `theta0`, `theta1` (los verderos parámetros)

```python
n = len(X)

# Escribe tu código para theta0_opt y theta1_opt
theta1_opt = 0
theta0_opt = 1

xrange = np.linspace(-1, 1.5, 100)
plt.scatter(X, y)
plt.plot(xrange, theta0_opt + theta1_opt * xrange, color="tab:orange")
```

## La Regresión Lineal (v.2)

¿Qué sucede si tenemos más de un feature (variable independiente) por cada observación?
Supongamos que tenemos $n$ observaciones, cada una con $m$ features.

x1|x2|...|xn|y
-|-|-|-|-
$x_{1,1}$|$x_{1,2}$|$\ldots$|$x_{1,n}$|$y_1$
$x_{2,1}$|$x_{2,2}$|$\ldots$|$x_{2,n}$|$y_2$
$x_{3,1}$|$x_{3,2}$|$\ldots$|$x_{3,n}$|$y_3$
$\vdots$|$\vdots$|$\ldots$|$\vdots$|$\vdots$
$x_{m,1}$|$x_{m,2}$|$\ldots$|$x_{m,n}$|$y_{m}$

Donde $x_{i,j} = x^{(i)}_j$ indica la $j$-ésima característica de la $i$-ésima observación.

En este caso, el problema sería encontrar $\theta_1, \theta_2, \ldots, \theta_m$. Nuestra predicción para la $i$-ésima observación $\hat y$ estaría dada por
$$
    \begin{align}
    \hat y_i &= \theta_0 + \theta_1 x_{1} +  \theta_2 x_2 + \ldots + \theta_m x_m\\
            &= \theta^T x^i
    \end{align}
$$

Donde
$$
    \theta = \begin{bmatrix}
                \theta_0 \\
                \theta_1 \\
                \theta_2 \\
                \vdots  \\
                \theta_m
             \end{bmatrix}
$$

$$
    x_i = \begin{bmatrix}
                1 \\
                x_{1,i} \\
                x_{2,i} \\
                \vdots  \\
                x_{m,i}
             \end{bmatrix}
$$

Definimos nuestro *training* dataset como la siguiente representación matricial:
$$
    X = \begin{bmatrix}
        – x_1 – \\
        – x_2 – \\
        \vdots \\
        – x_m – \\
        \end{bmatrix}
        = \begin{bmatrix}
                1 & x_{1,1} & x_{1,2} & \ldots & x_{1,m}\\
                1 & x_{2,1} & x_{2,2} & \ldots & x_{2,m}\\
                \vdots & \vdots & \vdots & \ddots & \vdots\\
                1 & x_{n,1} & x_{n,2} & \ldots & x_{n,m}\\
          \end{bmatrix}
$$

siendo los targets:
$$
    y = \begin{bmatrix} y_1 & y_2 & \ldots & y_n \end{bmatrix}
$$

Para un problema de **regresión lineal generalizada**, el vector $\theta^* = [\theta_0, \theta_1, \ldots, \theta_m]$ de valores que minimizan $J(\theta_0, \theta_1, \ldots, \theta_m)$ dado por

$$
    \theta^* = (X^T X)^{-1} X^T y
$$

Ejercicio

Considerando

```python
Xx = np.c_[np.ones(len(X)), X]
```

y el vector $y$, usando la fórmula de la regresión lineal generalizada, valida que las entradas de este vector sean iguales a las variables `theta0_opt`, `theta1_opt`.

**Nota 1**: La matrix `Xx` añade una columna de $1$s para cálcular el *bias* o factor de corrección

**Nota 2**: Para calcular la inversa de una matriz puedes usar la función `np.linalg.inv` desde `numpy`

### Regresión Linal usando `scikit-learn`

Durante este curso estaremos usando la librería `scikit-learn` para ahorrarnos los detalles en la implementación de los modelos. `scikit-learn` está escrita sobre numpy y es usada por compañías y universidades en su implementación y análisis de modelos de ML.

El primer paso para usar cualquier modelo dentro de la librería scikit learn es importarlo desde la *familia* a la cuál pertenece. Para nuestro modelo de regresión lineal, el modelo es llamado `LinearRegressio` y se encuentra dentro de la familia `linear_model`.

```python
from sklearn.linear_model import LinearRegression

# Inicializamos una instancia de la clase LinearRegression
# Hacer esto permite modificar hyperparámetros del modelo a aprender
linreg = LinearRegression()
# Pasamos la mariz de features y el vector de targets
linreg.fit(X.reshape(-1, 1), y)
linreg.intercept_, linreg.coef_
```

    (0.30375826386053706, array([1.55643992]))

### Interpretando los resultados de una regresión lineal

Una regresión lineal, aunque un modelo relativamente simple, nos da la oportunidad de poder interpretar los resultados obtenidos al problema sobre el cual estamos trabajando.

Consideremos el dataset `california_housing` con el precio de las casas en california. De manera ingenua, consideraremos todas las variables dentro de la base de datos *as is* (recordemos que antes de entrenar cualquier modelo de machine learning, es importante analizar los datos para poder llegar a una conclusión más razonable sobre los resultados del modelo o como mejorar su desempeño)

```python
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
housing = fetch_california_housing()
print(housing.DESCR)
print(housing.feature_names)
```

    .. _california_housing_dataset:

    California Housing dataset
    --------------------------

    **Data Set Characteristics:**

        :Number of Instances: 20640

        :Number of Attributes: 8 numeric, predictive attributes and the target

        :Attribute Information:
            - MedInc        median income in block
            - HouseAge      median house age in block
            - AveRooms      average number of rooms
            - AveBedrms     average number of bedrooms
            - Population    block population
            - AveOccup      average house occupancy
            - Latitude      house block latitude
            - Longitude     house block longitude

        :Missing Attribute Values: None

    This dataset was obtained from the StatLib repository.
    http://lib.stat.cmu.edu/datasets/

    The target variable is the median house value for California districts.

    This dataset was derived from the 1990 U.S. census, using one row per census
    block group. A block group is the smallest geographical unit for which the U.S.
    Census Bureau publishes sample data (a block group typically has a population
    of 600 to 3,000 people).

    It can be downloaded/loaded using the
    :func:`sklearn.datasets.fetch_california_housing` function.

    .. topic:: References

        - Pace, R. Kelley and Ronald Barry, Sparse Spatial Autoregressions,
          Statistics and Probability Letters, 33 (1997) 291-297

    ['MedInc', 'HouseAge', 'AveRooms', 'AveBedrms', 'Population', 'AveOccup', 'Latitude', 'Longitude']

```python
housing.data
```

    array([[   8.3252    ,   41.        ,    6.98412698, ...,    2.55555556,
              37.88      , -122.23      ],
           [   8.3014    ,   21.        ,    6.23813708, ...,    2.10984183,
              37.86      , -122.22      ],
           [   7.2574    ,   52.        ,    8.28813559, ...,    2.80225989,
              37.85      , -122.24      ],
           ...,
           [   1.7       ,   17.        ,    5.20554273, ...,    2.3256351 ,
              39.43      , -121.22      ],
           [   1.8672    ,   18.        ,    5.32951289, ...,    2.12320917,
              39.43      , -121.32      ],
           [   2.3886    ,   16.        ,    5.25471698, ...,    2.61698113,
              39.37      , -121.24      ]])

```python
housing.target
```

    array([4.526, 3.585, 3.521, ..., 0.923, 0.847, 0.894])

```python
linmodel = LinearRegression()
xtrain, xtest, ytrain, ytest = train_test_split(housing.data, housing.target,
                                                test_size=0.2, random_state=314)
linmodel.fit(xtrain, ytrain)
print(f"R2 del modelo: {linmodel.score(xtrain, ytrain):0.5f}")
```

    R2 del modelo: 0.60665

Ejercicio

Dados los resultados de nuestra regresión lineal sobre los datos `housing`, interpreta los resultados.

```python
for coef, lab in zip(linmodel.coef_, housing.feature_names):
    print(f"{lab:>{len(max(housing.feature_names))}}: {coef * 100_000:,.2f}")
```

        MedInc: 43,619.53
      HouseAge: 980.06
      AveRooms: -10,515.12
     AveBedrms: 61,910.56
    Population: -0.35
      AveOccup: -408.19
      Latitude: -41,988.09
     Longitude: -43,344.95

### Generalización: *overfitting* v.s. *underfitting*

```python
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.metrics import mean_squared_error

degrees = np.arange(1, 8)
seed(31415)
nobs = 50

Xtrain = np.linspace(0, 1, nobs) + randn(nobs) * 0.15
ytrain = 0.3 + 0.7 * Xtrain - 2.3 * Xtrain ** 2 - 1.4 * Xtrain ** 3 + randn(nobs) * 0.35

Xtest = np.linspace(0, 1, nobs) + randn(nobs) * 0.15
ytest = 0.3 + 0.7 * Xtest - 2.3 * Xtest ** 2 - 1.4 * Xtest ** 3 + randn(nobs) * 0.35

Xtrain = Xtrain.reshape(-1, 1)
Xtest = Xtest.reshape(-1, 1)

train_cost = []
test_cost = []
for deg in degrees:
    #sklearn pipeline: transformamos la variable x en
    # un polinomio de grado deg, seguido de una regresión lineal
    lreg_pipe = Pipeline([
        ("poly", PolynomialFeatures(deg)),
         ("lreg", LinearRegression())
    ])

    lreg_pipe.fit(Xtrain, ytrain)
    train_preds = lreg_pipe.predict(Xtrain)
    test_preds = lreg_pipe.predict(Xtest)

    train_cost.append(mean_squared_error(ytrain, train_preds))
    test_cost.append(mean_squared_error(ytest, test_preds))
```

Hacer un modelo para mejorar el costo sobre el **training** dataset no necesariamente implica una mejor generalización a datos aún no observados.

La habilidad de un modelo para aprender una función compleja de los datos se llama **capacidad**. A mayor capacidad, mayor complejidad de los datos podemos capturar. Mucha capacidad y el modelo sufre de *overfitting*, muy poca y el modelo sufre de *underfitting* En general:

    <bold> Complejidad $\not\Rightarrow$ Precisión </bold>

Supongamos que
$$
    y = 0.3 + 0.7 x - 2.3 x^2 - 1.4 x ^3 + \varepsilon
$$

```python
def plot_costs(degrees, train, test):
    plt.plot(degrees, train)
    plt.scatter(degrees, train, label="train")
    plt.plot(degrees, test)
    plt.scatter(degrees, test, label="test")
    plt.grid(alpha=0.5)
    plt.title("Costo del modelo como función de un polinomío")
    plt.xlabel("Degree")
    plt.ylabel(r"$J(\theta)$")
    plt.legend()
```

```python
plot_costs(degrees, train_cost, test_cost)
```

max_stars_repo_name: dlcruzsergio/predict-home-value

max_issues_repo_name: dlcruzsergio/predict-home-value

max_forks_repo_name: dlcruzsergio/predict-home-value


text:
# numpy

* faster computing than native python (C-wrapper), **30x faster** `np.array` manipulation than `list`
* favourite in data processing tasks, AI (keras)
* image (via `cv2`) as numpy 2D array
* pandas for data analysis

```python
import numpy as np
```

```python
# if it's not installed yet, let's use "magic command" `!`
!pip install numpy
```

    Collecting numpy
      Using cached numpy-1.20.1-cp38-cp38-win_amd64.whl (13.7 MB)
    Installing collected packages: numpy
    Successfully installed numpy-1.20.1

    WARNING: You are using pip version 20.1.1; however, version 21.0.1 is available.
    You should consider upgrading via the 'c:\prace\python-workshop-org\venv\scripts\python.exe -m pip install --upgrade pip' command.

```python
numbers = range(5)
a = np.array(numbers)
print(a)

b = np.arange(5)
print(b)
```

    [0 1 2 3 4]
    [0 1 2 3 4]

```python
# compare each item individually
print(a == b)
```

    [ True  True  True  True  True]

```python
# compare list content as singe object
list([1,2,3,4,5]) == list([1,2,3,4,5])
```

    True

# range vs linspace

```python
numbers = list(range(0, 10))
print(numbers)
```

    [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]

```python
numbers = np.linspace(0, 10)
print(numbers)
```

    [ 0.          0.20408163  0.40816327  0.6122449   0.81632653  1.02040816
      1.2244898   1.42857143  1.63265306  1.83673469  2.04081633  2.24489796
      2.44897959  2.65306122  2.85714286  3.06122449  3.26530612  3.46938776
      3.67346939  3.87755102  4.08163265  4.28571429  4.48979592  4.69387755
      4.89795918  5.10204082  5.30612245  5.51020408  5.71428571  5.91836735
      6.12244898  6.32653061  6.53061224  6.73469388  6.93877551  7.14285714
      7.34693878  7.55102041  7.75510204  7.95918367  8.16326531  8.36734694
      8.57142857  8.7755102   8.97959184  9.18367347  9.3877551   9.59183673
      9.79591837 10.        ]

```python
numbers = np.linspace(0, 10, 50)
print(numbers)
```

    [ 0.          0.20408163  0.40816327  0.6122449   0.81632653  1.02040816
      1.2244898   1.42857143  1.63265306  1.83673469  2.04081633  2.24489796
      2.44897959  2.65306122  2.85714286  3.06122449  3.26530612  3.46938776
      3.67346939  3.87755102  4.08163265  4.28571429  4.48979592  4.69387755
      4.89795918  5.10204082  5.30612245  5.51020408  5.71428571  5.91836735
      6.12244898  6.32653061  6.53061224  6.73469388  6.93877551  7.14285714
      7.34693878  7.55102041  7.75510204  7.95918367  8.16326531  8.36734694
      8.57142857  8.7755102   8.97959184  9.18367347  9.3877551   9.59183673
      9.79591837 10.        ]

```python
np.linspace?
```

# matrix shape

```python
matrix = np.arange(60)
matrix
```

    array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
           17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
           34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50,
           51, 52, 53, 54, 55, 56, 57, 58, 59])

```python
matrix.shape
```

    (60,)

```python
matrix.reshape(4, 15)
```

    array([[ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14],
           [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29],
           [30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44],
           [45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59]])

```python
matrix.reshape(6, 10)
```

    array([[ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9],
           [10, 11, 12, 13, 14, 15, 16, 17, 18, 19],
           [20, 21, 22, 23, 24, 25, 26, 27, 28, 29],
           [30, 31, 32, 33, 34, 35, 36, 37, 38, 39],
           [40, 41, 42, 43, 44, 45, 46, 47, 48, 49],
           [50, 51, 52, 53, 54, 55, 56, 57, 58, 59]])

```python
# original shape is just 1D vector
matrix.reshape(60,)
```

    array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
           17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
           34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50,
           51, 52, 53, 54, 55, 56, 57, 58, 59])

```python
m = matrix.reshape(6, 10)
print(m)
print(m.shape)
```

    [[ 0  1  2  3  4  5  6  7  8  9]
     [10 11 12 13 14 15 16 17 18 19]
     [20 21 22 23 24 25 26 27 28 29]
     [30 31 32 33 34 35 36 37 38 39]
     [40 41 42 43 44 45 46 47 48 49]
     [50 51 52 53 54 55 56 57 58 59]]
    (6, 10)

#### @ show trick with double click to left header of cell

```python
# transform to 3D
print(matrix.reshape(6, 5, 2))
```

    [[[ 0  1]
      [ 2  3]
      [ 4  5]
      [ 6  7]
      [ 8  9]]

     [[10 11]
      [12 13]
      [14 15]
      [16 17]
      [18 19]]

     [[20 21]
      [22 23]
      [24 25]
      [26 27]
      [28 29]]

     [[30 31]
      [32 33]
      [34 35]
      [36 37]
      [38 39]]

     [[40 41]
      [42 43]
      [44 45]
      [46 47]
      [48 49]]

     [[50 51]
      [52 53]
      [54 55]
      [56 57]
      [58 59]]]

# slicing

```python
print(m)
```

    [[ 0  1  2  3  4  5  6  7  8  9]
     [10 11 12 13 14 15 16 17 18 19]
     [20 21 22 23 24 25 26 27 28 29]
     [30 31 32 33 34 35 36 37 38 39]
     [40 41 42 43 44 45 46 47 48 49]
     [50 51 52 53 54 55 56 57 58 59]]

```python
# slice by 1. dimension (row)
m[1]
```

    array([10, 11, 12, 13, 14, 15, 16, 17, 18, 19])

```python
m[1:4]
```

    array([[10, 11, 12, 13, 14, 15, 16, 17, 18, 19],
           [20, 21, 22, 23, 24, 25, 26, 27, 28, 29],
           [30, 31, 32, 33, 34, 35, 36, 37, 38, 39]])

```python
# slice by 2. dimension (column)
m[:, 2]
```

    array([ 2, 12, 22, 32, 42, 52])

```python
# subset (indexing from 0!)
m[2:4, 3:7]
```

    array([[23, 24, 25, 26],
           [33, 34, 35, 36]])

```python
# subset (indexing from 0!)
m[2:4, 3:7]
```

    array([[23, 24, 25, 26],
           [33, 34, 35, 36]])

#### @ slice by axes with comma, or [ ]

```python
m[2][5]
```

    25

```python
m[2,5]
```

    25

# modify matrix by assignment of single value
#### @ single value replacement

```python
m[5,2] = 999
print(m)
```

    [[  0   1   2   3   4   5   6   7   8   9]
     [ 10  11  12  13  14  15  16  17  18  19]
     [ 20  21  22  23  24  25  26  27  28  29]
     [ 30  31  32  33  34   0   0  37  38  39]
     [ 40  41  42  43  44   0   0  47  48  49]
     [ 50  51 999  53  54  55  56  57  58  59]]

```python
number = m[5,2]
print(number)
```

    999

```python
number = 123
print(number)
```

    123

```python
# and it's not modified (as expected)
print(m)
```

    [[  0   1   2   3   4   5   6   7   8   9]
     [ 10  11  12  13  14  15  16  17  18  19]
     [ 20  21  22  23  24  25  26  27  28  29]
     [ 30  31  32  33  34   0   0  37  38  39]
     [ 40  41  42  43  44   0   0  47  48  49]
     [ 50  51 999  53  54  55  56  57  58  59]]

# modify matrix by block assignment

```python
# maybe there is prettier way
m = np.arange(60).reshape(6, 10)
print(m)
```

    [[ 0  1  2  3  4  5  6  7  8  9]
     [10 11 12 13 14 15 16 17 18 19]
     [20 21 22 23 24 25 26 27 28 29]
     [30 31 32 33 34 35 36 37 38 39]
     [40 41 42 43 44 45 46 47 48 49]
     [50 51 52 53 54 55 56 57 58 59]]

### np.zeros

```python
zeros = np.zeros((2,2))
print(zeros)
```

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

```python
print(m[3:5, 5:7])
```

    [[35 36]
     [45 46]]

note: same shape

```python
m[3:5, 5:7] = zeros
print(m)
```

    [[ 0  1  2  3  4  5  6  7  8  9]
     [10 11 12 13 14 15 16 17 18 19]
     [20 21 22 23 24 25 26 27 28 29]
     [30 31 32 33 34  0  0 37 38 39]
     [40 41 42 43 44  0  0 47 48 49]
     [50 51 52 53 54 55 56 57 58 59]]

# matrix views - replacement within submatrix

* slicing returns *views* (=shallow copy)
* `.copy()` to get deep copy

```python
m = np.arange(60).reshape(6, 10)
print(m)
```

    [[ 0  1  2  3  4  5  6  7  8  9]
     [10 11 12 13 14 15 16 17 18 19]
     [20 21 22 23 24 25 26 27 28 29]
     [30 31 32 33 34 35 36 37 38 39]
     [40 41 42 43 44 45 46 47 48 49]
     [50 51 52 53 54 55 56 57 58 59]]

```python
submatrix = m[0:3, 0:3]
print(submatrix)
```

    [[ 0  1  2]
     [10 11 12]
     [20 21 22]]

```python
submatrix[1,1] = 9999
print(submatrix)
```

    [[   0    1    2]
     [  10 9999   12]
     [  20   21   22]]

```python
print(m)
```

    [[   0    1    2    3    4    5    6    7    8    9]
     [  10 9999   12   13   14   15   16   17   18   19]
     [  20   21   22   23   24   25   26   27   28   29]
     [  30   31   32   33   34   35   36   37   38   39]
     [  40   41   42   43   44   45   46   47   48   49]
     [  50   51   52   53   54   55   56   57   58   59]]

# Matrix basics

### broadcasting

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

    [[1 2]
     [3 4]]

```python
a + 100
```

    array([[101, 102],
           [103, 104]])

```python
a * 5
```

    array([[ 5, 10],
           [15, 20]])

#### `*` will still apply broadcasting!

```python
# np.array([1,5]) works as well
# a * np.array([1, 5])
(a * np.array([1, 5]))
```

    array([[ 1, 10],
           [ 3, 20]])

# matrix multiplication

```python
a.dot([1,5])
```

    array([11, 23])

matrix muliplication has special operator in python `@`

```python
a @ np.array([1, 5])
```

    array([11, 23])

## transposing matrix

```python
a
```

    array([[1, 2],
           [3, 4]])

```python
a.T
```

    array([[1, 3],
           [2, 4]])

# compare matrixes

```python
a = np.arange(5)
b = np.arange(5)

print(f'{a=}')
print(f'{b=}')
```

    a=array([0, 1, 2, 3, 4])
    b=array([0, 1, 2, 3, 4])

#### @ compare as single entity

```python
list(a) == list(b)
```

    True

#### @ item by item comparison

```python
a == b
```

    array([ True,  True,  True,  True,  True])

#### @ comapre operators as 'mask'

```python
a >= 2
```

    array([False, False,  True,  True,  True])

```python
a[a >= 2]
```

    array([2, 3, 4])

#### @ is same as mask from list of bools

```python
mask = [False, False, True, True, True]
a[mask]
```

    array([2, 3, 4])

# `.flatten` vs `.ravel`

* `.flatten` creates deep copy
* `.ravel` work as *view* (more memory efficient)

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

    [[1 2]
     [3 4]]

```python
changed = matrix.ravel()
print(changed)
```

    [1 2 3 4]

```python
changed = matrix.flatten()
print(changed)
```

    [1 2 3 4]

```python
changed[3] = 999
print(changed)
print("")
print(matrix)
```

    [  1   2   3 999]

    [[1 2]
     [3 4]]

# save / load `np.ndarray`

* `np.save`, `np.load` - save dimension, datatype, ...
* `np.savetxt`, `np.loadtxt` - save as CSV

```python
vector = np.arange(10).reshape(2,5)
print(vector)
np.save('export.npy', vector)
np.savetxt('export.csv', vector)
```

    [[0 1 2 3 4]
     [5 6 7 8 9]]

```python
np.load('export.npy')
```

    array([[0, 1, 2, 3, 4],
           [5, 6, 7, 8, 9]])

```python
np.loadtxt('export.csv')
```

    array([[0., 1., 2., 3., 4.],
           [5., 6., 7., 8., 9.]])

```python
%load export.csv
```

```python
m = matrix.reshape(6, 10)
m
```

    array([[ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9],
           [10, 11, 12, 13, 14, 15, 16, 17, 18, 19],
           [20, 21, 22, 23, 24, 25, 26, 27, 28, 29],
           [30, 31, 32, 33, 34, 35, 36, 37, 38, 39],
           [40, 41, 42, 43, 44, 45, 46, 47, 48, 49],
           [50, 51, 52, 53, 54, 55, 56, 57, 58, 59]])

# let numpy to solve linear equations

It's an individual task as homework 😉

# sympy
https://docs.sympy.org/latest/index.html

* symbolic computation
  * simplyfing expressions
  * variable substitions
* integrals
* derivatives
* LaTeX output support (`sympy.init_printing(use_latex=False, use_unicode=False)`)

```python
import sympy
x = sympy.symbols('x')
```

```python
expr = 2*x + 3
expr
```

$\displaystyle 2 x + 3$

```python
expr.subs(x, 2)
```

$\displaystyle 7$

```python
a, b, c = sympy.symbols('a b c')
expr = (x + 1)*(x-1)
expr
```

$\displaystyle \left(x - 1\right) \left(x + 1\right)$

```python
sympy.simplify(expr)
```

$\displaystyle x^{2} - 1$

```python
sympy.Integral(1/x) + sympy.sqrt(x + 1)
```

max_stars_repo_name: mpavlase/python-workshop

max_issues_repo_name: mpavlase/python-workshop

max_forks_repo_name: mpavlase/python-workshop