ryamadaの遺伝学・遺伝統計学メモ このページをアンテナに追加 RSSフィード

数学・コンピュータ関連の姉妹ブログ『ryamadaのコンピュータ・数学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

2017-01-22 統計遺伝学分野で生き抜くための数学〜精選8〜

[][][][]期待値計算のための微積分

---
title: "Calculus1 Expected value and basics of differentiation and integration"
author: "ryamada"
date: "2017年1月22日"
output: 
  html_document:
    toc: true
    toc_depth: 6
    number_section: true
---

```{r setup, include=FALSE}
library(knitr)
opts_chunk$set(echo = TRUE)
library(rgl)
knit_hooks$set(rgl = hook_webgl)
```

# 平均 Average

$$
m = \frac{1}{n} \sum_{i=1}^n x_i
$$

```{r}
n <- 10
x <- sample(1:3,n,replace=TRUE)
x
mean(x)
1/n * sum(x)
```

# 重み付き平均 Weighted average

$$
m_w = \sum_{j=1}^k v_j \times Pr(j)
$$
```{r}
x
tabulate(x)
w <- tabulate(x)/n
w
v <- sort(unique(x))
v
mw <- sum(v * w)
mw
```

# 期待値 Expected value

## サイコロの目の数の期待値 Expected value of dice

```{r}
v <- 1:6
p <- rep(1/6,6)
sum(v*p)
```

出来の悪いサイコロの期待値 Dice in bad condition

```{r}
v <- 1:6
p <- rep(1/6,6) + rnorm(6) * 0.1
p <- p/sum(p)
plot(v,p,type="h")
p
sum(p)
sum(v*p)
```

## 二項分布の期待値は $np$ Expected value of binomial distribution : $np$

$$
(p + (1-p)) ^n = 1^n = 1 = \sum_{i=0}^n \begin{pmatrix}n\\i \end{pmatrix} p^i (1-p)^i
$$
```{r}
p0 <- 0.3
n <- 10
i <- 0:n
i.inv <- i[(n+1):1]
choose(n,i)
p <- choose(n,i) * p0^i * (1-p0)^i.inv
plot(i,p,type="h")
sum(p)
sum(p*i)

n * p0
```

## ベータ分布の正規化 Normalization of beta distribution

成功・失敗が、n回とm回だったとき、成功率がpである尤度は n successes and m failures. Likelihood of success rate p is proportional to ;

$$
p^n (1-p)^m
$$
に比例する。

With $h(n,m)$ below,
$$
h(n,m)=\int_0^1 p^n(1-p)^m dp
$$

とおけば、the following equaion follows;

$$
\int_0^1 \frac{1}{h(n,m)}p^n (1-p)^m dp =1
$$
となるから

The following is the likelihood function of success rate p when n successes and m failures.
$$
\frac{1}{h(n,m)}p^n (1-p)^m
$$

が成功n回、失敗m回のときの成功率pの尤度関数。

$p^n(1-p)^m$が関数の形を決め、$h(n,m)$は積分が1となるように正規化しているので、$h(n,m)$によって(尤度)関数を正規化する、と言う。

$p^n (1-p)^m$ determines its shape and $h(n,m)$ normalizes its integration to 1.


## $h(n,m)=\int_0^1 p^n(1-p)^m dp$を計算してみる Calculation of $h(n,m)$

### n=1, m=0
$$
h(1,0) = \int_0^1 p dp
$$

```{r}
p <- seq(from=0,to=1,length=100)
h <- p
plot(p,h,type="l")
```

$h(1,0)$は面積として計算できる。 Area of $h(1,0)$ is given geometrically. 

$$
\frac{1}{2} \times 1 \times 1 = \frac{1}{2}
$$


積分するなら Integration;

$$
\frac{d}{dx} x^2 = 2x\\
\frac{d}{dx} \frac{1}{2} x^2 = x
$$
を使って、

$$
x^2 + C = \int 2x dx\\
\frac{1}{2} (x^2+C) = \int x dx
$$

から、

$$
\int_0^1 p dp = [\frac{1}{2} x^2]^1_0 = \frac{1}{2} (1^2-0^2)=\frac{1}{2}
$$
となる。

結局、n=1,m=0のときの尤度関数は The likelihood function when n=1 and m=0;
$$
\frac{1}{h(1,0)}p^1 (1-p)^0 = \frac{1}{\frac{1}{2}}p = 2p
$$

# Exercise 1

## Exercise 1-1
n=1, m=1の場合、二項観察の尤度関数は Likelihood function for binomial observation n=1 and m=1 
$$
\frac{1}{h(1,1)} p(1-p) = \frac{1}{h(1,1)}(p-p^2)
$$
$$
h(1,1) = \int_0^1 p-p^2 dp
$$

を求めたい。

$f(x) = x^1 (1-x)^1$ のグラフを描け。Draw $f(x) = x^1(1-x)^1$.

## Exercise 1-2

[0,1]区間を、k等分してその小区間ごとの面積を近似的に計算し、その和を[0,1]の範囲の$p-p^2$の面積とみなすこととする。
第i小区間の面積を、長方形の面積とみなして、計算し、kを、1,2,...,100と変化させ、その様子をプロットせよ。ただし、長方形は幅$\frac{1}{k}$、高さはその小区間の両端の$p-p^2$の値の平均値とせよ。

Divide the interval [0,1] into k evenly. Calculate subintervals' area approximately and sum them which is approcimation of the area under the curve.
The area of the i-th subinterval should be considered a rectangle whose width is $/frac{1}{k}$ and its hight is the average of the hights of the both ends of the rectangle. Calculate and plot for k=1,2,...,100.

## Exercise 1-3

$\frac{d}{dx}x^2 = 2x, \frac{d}{dx}x^3=3x^2$を使って$h(1,1)$を求め、近似で求めた値と比較せよ。

Integrate the function and compare the value with the approximation above.

## Exercise 1-4

n=2,m=3の場合の$p^n(1-p)^3$を展開し、n=1,m=1の場合と同様のことをせよ

Do the same for n=2 and m=3.

## Exercise 1-5

指数分布の期待値は $\frac{1}{\lambda}$であると言う。このことを、離散的な計算をすることで確認せよ。

The expected value of exponential distribution is $\frac{1}{\lambda}$. Calculate its expected value discretely.

$$
Pr(x) = \lambda \times e^{-\lambda x}
$$

## Exercise 1-6
微分積分の基礎技術 Basic skills of calculus 

Go through the every item in the page 
https://en.wikibooks.org/wiki/Calculus/Differentiation/Basics_of_Differentiation/Solutions#Find_The_Derivative_By_Definition 

2016-12-31 ぱらぱらめくる『量子力学で生命の謎を解く』

[][]ぱらぱらめくる『量子力学で生命の謎を解く』

  • 目次
    • 1 はしがき
    • 2 生命とは何か?
    • 3 生命のエンジン
    • 4 量子のうなり
    • 5 ニモの家を探せ
    • 6 チョウ、ショウジョウバエ、量子のコマドリ
    • 7 量子の遺伝子
    • 8 心
    • 9 生命の起源
    • 10 量子生物学--嵐の淵の生命
    • エピローグ 量子革命
  • 1 はしがき
    • 量子力学とは何か
    • 粒子・波動の2面性
    • 量子トンネル効果
    • 量子重ね合わせ
    • 確率的状態
    • 量子もつれ
    • ミクロな世界だがマクロにも影響。生物を含む日常生活でも活用されている(MRIの撮像は量子力学反応を生体にて観測すること。トリの方向磁気感覚も磁場に敏感な分子反応基の量子化学的状態と関連か)
  • 2 生命とは何か?
    • 無秩序から秩序へ
    • 熱力学法則に従わない
    • 熱力学法則は、「たくさんの分子」x「それなりの広さ」x「自由な運動」→「平均が重要」という考え方
    • それに従わないときが無秩序から秩序を作り出しそう
      • 「少ない分子」、「狭い空間」、「自由に運動できない」→確率的現象
  • 3 生命のエンジン
    • 酵素の話
    • 触媒とは何か
    • 反応閾値(平均としての)を下げること〜古典熱力学的説明
    • 量子トンネル効果による、低温での反応進行〜量子力学的説明
  • 4 量子のうなり
    • 植物は量子コンピュータである
    • 量子コンピュータは、「観測するまでは、色々な状態を同時に取っている」という量子的現象を使って計算するコンピュータ
    • 「観測しない」ということは、「分子・量子」に「何も触れない・反応を仕掛けない」という状態
    • そんな状態は、超低温や、超低密度でしか起きないのが普通だが、「植物体」はそんな環境を常温で実現しているのでは、と言う話
    • 光合成において、光子が、光合成分子の反応中心に辿り着けるのは、たくさんある通り道の「最適解」を量子コンピュータ的な探索をすることで見つけ出すから
  • 5 ニモの家を探せ
    • 嗅覚
    • 嗅覚受容体はそこそこの種類の遺伝子を組み合わせることで膨大な化学分子の認識を可能にしている
    • 鍵と鍵穴的な対応関係〜形状での分子認識
    • 化学分子に光を当てて、その振動スペクトルを検知する〜振動による分子認識
    • 嗅覚受容体は、両方を併せて使う
  • 6 チョウ、ショウジョウバエ、量子のコマドリ
    • 磁気コンパス
    • 量子もつれは化学反応に影響を与える
    • トリの磁気感覚では量子もつれ状態にある遊離基のペアが関わる
  • 7 量子の遺伝学
    • 遺伝子は「非周期的結晶」
    • DNA二重らせんの塩基対は量子効果
    • 突然変異は、量子的化学反応でも起きうる
    • 遺伝子が離散情報分子であること自身に量子性を求めている章?
  • 8 心
    • 意識が、最終的に個体の行動を起こす、その様子を、脳の量子コンピュータ性で説明しようとする
  • 9 生命の起源
    • 自己複製という反応が始まったとき
  • 10 漁師生物学--嵐の縁の生命
    • 自己組織化反応は「縁」で起きている
  • エピローグ 量子革命

2016-12-28 統計遺伝学分野で生き抜くための数学〜精選7.5〜

[][][][][]最小二乗法のための微分

---
title: "Differentiate Matrix Formula 線形代数式を微分する"
author: "ryamada"
date: "2016年12月27日"
output: 
  html_document:
    toc: true
    toc_depth: 6
    number_section: true
---

```{r setup, include=FALSE}
library(knitr)
opts_chunk$set(echo = TRUE)
library(rgl)
knit_hooks$set(rgl = hook_webgl)
```

# 最小二乗法

$$
y \sim X\mathbf{a}
$$
なる関係式があり、$X$の行数が列数より多いとき、

$$
\sum_{i=1}^n (y_i-\hat{y_i})^2 = (y-X\mathbf{a})^T \cdot (y-X\mathbf{a})
$$
を最小にするような$\hat{\mathbf{a}}$が、推定値であった。

そして
$$
\hat{\mathbf{a}} = (X^TX)^{-1}X^Ty
$$

で求まるのであった。

# 偏微分方程式

$$
f(\mathbf{a})=f((a_1,...,a_m))=\sum_{i=1}^n (y_i-\hat{y_i})^2 = (y-X\mathbf{a})^T \cdot (y-X\mathbf{a})
$$
は、$y$,$X$が与えられているとき、$m$個の変数$(a_1,...,a_m)$の二次式となっているスカラー関数である。

今、この$f(\mathbf{a})$を最小にする$\mathbf{a}=(a_1,...,a_m)$とは、
$$
\frac{\partial f(\mathbf{a})}{\partial a_i}=0
$$
がすべての$a_i$について成り立つ場合である。

スカラーを返す$z^Tw$は$z^Tw = w z^T$であることを使うと

$$
f(\mathbf{a}) = (y-X\mathbf{a})^T \cdot (y-X\mathbf{a})
 = y^Ty -2(X\mathbf{a})^Ty+(X\mathbf{a})^T \cdot (X\mathbf{a})
$$
さらに変形して

$$
f(\mathbf{a}) = y^Ty -2(X^Ty)\cdot \mathbf{a} + \mathbf{a}^TX^TX\mathbf{a}
$$

これを$\mathbf{a}$で偏微分する。

第1項は$0$

第2項は、各成分に$-2y^TX$の各成分が残る。

第3項は、$X^TX=Z$と置くと

$$
\mathbf{a}^T Z \mathbf{a} = \sum_{i=1}^m \sum_{j=1}^m a_i a_j Z_{ij}
$$

であり、その$a_k$による偏微分は

$$
\frac{\partial}{\partial a_k} (\sum_{i=1}^m \sum_{j=1}^m a_i a_j Z_{ij}) = 2\sum_{i=1}^m Z_{ik}a_i
$$

となり、$k=1,...,m$について合わせると
$$
2Z \mathbf{a}
$$
となるから、
結局、

$$
0 - 2X^Ty + 2Z\mathbf{a} = 2(X^TX) \mathbf{a} - 2X^Ty = \mathbf{0}
$$

となり、

$$
(X^TX)\mathbf{a} = X^T y
$$
が得られる。

ここから、$X^TX$に逆行列があるときは

$$
\mathbf{a} = (X^TX)^{-1}X^T y
$$

が得られる。

[][][][][]最小二乗法のための微分

---
title: "Differentiate Matrix Formula 線形代数式を微分する"
author: "ryamada"
date: "2016年12月27日"
output: 
  html_document:
    toc: true
    toc_depth: 6
    number_section: true
---

```{r setup, include=FALSE}
library(knitr)
opts_chunk$set(echo = TRUE)
library(rgl)
knit_hooks$set(rgl = hook_webgl)
```

# 最小二乗法

$$
y \sim X\mathbf{a}
$$
なる関係式があり、$X$の行数が列数より多いとき、

$$
\sum_{i=1}^n (y_i-\hat{y_i})^2 = (y-X\mathbf{a})^T \cdot (y-X\mathbf{a})
$$
を最小にするような$\hat{\mathbf{a}}$が、推定値であった。

そして
$$
\hat{\mathbf{a}} = (X^TX)^{-1}X^Ty
$$

で求まるのであった。

# 偏微分方程式

$$
f(\mathbf{a})=f((a_1,...,a_m))=\sum_{i=1}^n (y_i-\hat{y_i})^2 = (y-X\mathbf{a})^T \cdot (y-X\mathbf{a})
$$
は、$y$,$X$が与えられているとき、$m$個の変数$(a_1,...,a_m)$の二次式となっているスカラー関数である。

今、この$f(\mathbf{a})$を最小にする$\mathbf{a}=(a_1,...,a_m)$とは、
$$
\frac{\partial f(\mathbf{a})}{\partial a_i}=0
$$
がすべての$a_i$について成り立つ場合である。

スカラーを返す$z^Tw$は$z^Tw = w z^T$であることを使うと

$$
f(\mathbf{a}) = (y-X\mathbf{a})^T \cdot (y-X\mathbf{a})
 = y^Ty -2(X\mathbf{a})^Ty+(X\mathbf{a})^T \cdot (X\mathbf{a})
$$
さらに変形して

$$
f(\mathbf{a}) = y^Ty -2(X^Ty)\cdot \mathbf{a} + \mathbf{a}^TX^TX\mathbf{a}
$$

これを$\mathbf{a}$で偏微分する。

第1項は$0$

第2項は、各成分に$-2y^TX$の各成分が残る。

第3項は、$X^TX=Z$と置くと

$$
\mathbf{a}^T Z \mathbf{a} = \sum_{i=1}^m \sum_{j=1}^m a_i a_j Z_{ij}
$$

であり、その$a_k$による偏微分は

$$
\frac{\partial}{\partial a_k} (\sum_{i=1}^m \sum_{j=1}^m a_i a_j Z_{ij}) = 2\sum_{i=1}^m Z_{ik}a_i
$$

となり、$k=1,...,m$について合わせると
$$
2Z \mathbf{a}
$$
となるから、
結局、

$$
0 - 2X^Ty + 2Z\mathbf{a} = 2(X^TX) \mathbf{a} - 2X^Ty = \mathbf{0}
$$

となり、

$$
(X^TX)\mathbf{a} = X^T y
$$
が得られる。

ここから、$X^TX$に逆行列があるときは

$$
\mathbf{a} = (X^TX)^{-1}X^T y
$$

が得られる。

2016-12-27 統計遺伝学分野で生き抜くための数学〜精選7〜

[][][][][]行列で最善の解

---
title: "Best answer with Matrix: Moore-Penrose Pseudo-inverse 行列で最善の解: ムーアペンローズ疑似逆行列"
author: "ryamada"
date: "2016年12月27日"
output: 
  html_document:
    toc: true
    toc_depth: 6
    number_section: true
---

```{r setup, include=FALSE}
library(knitr)
opts_chunk$set(echo = TRUE)
library(rgl)
knit_hooks$set(rgl = hook_webgl)
```

# 連立方程式の解と最小二乗法による線形回帰

変数の数と等式の数が一致しているとき、うまく逆行列が取れれば
$$
y = X \mathbf{a}
$$
はきっちりと解けて$\mathbf{a}$が一意に求まる。

その際、
$$
\mathbf{a} = X^{-1} y = (X^TX)^{-1}X^T y
$$
であった。

線形回帰の場合は、サンプルの数(ベクトル$y$の長さ:$n$とする)に対して、説明変数の数$X$の列数($m$とする)が小さいので

$$
y \sim X \mathbf{a}
$$
として、
$$
\sum_{i=1}^n (y_i-\hat{y_i})^2
$$
が最小になるような$\mathbf{a}$を推定するわけだが、そのとき
$$
\mathbf{a} = (X^T X)^{-1} X^{-1}
$$
で求まるのだった。

結局、行列$X$の列数によらず、$\mathbf{a}$について、「これが一番」という答えが
$$
\mathbf{a} = (X^TX)^{-1}X^T y
$$
で求まる。

ここで、「これが一番」というのは、最小二乗が0であったり、「かっちり連立方程式が解ける」ということだったりした。

# 連立方程式でも線形回帰でもない場合

行列$X$の行数$n$,列数$m$について

+ $n=m$の場合:連立方程式を解く

+ $n>m$の場合:線形回帰をする

だった。

網羅されていないのは

+ $n<m$の場合である


この場合にも、うまく行く方法があり、これをムーアペンローズ疑似逆行列と読んだり、一般化逆行列と読んだり、疑似逆行列と呼んだりする。


$n < m$の場合には、$(X^TX)^{-1}$がうまく計算されないので、別の式を用いる。

Xの行数が列数より小さい場合には、転置して考える。
転置($X'=X^T$)すると、そのムーアペンローズ疑似逆行列は、$(X'^TX')^{-1}X'^T$と計算できるので、それを再度転置すれとうまく行く($((X'^TX')^{-1}X'^T)^T$)。

結局、覚えるべきは

+ $n \ge m$ のときは $gen.inv(X)=(X^TX)^{-1}X^T$

+ $n < m$のときは、$X_{gen.inv} = (gen.inv(X^T))^T$

$n < m$の場合を$X^T(X^TX)^{-1}$と書くこともある。


```{r}
n <- 3
m <- 2 # m < n
X <- matrix(rnorm(n*m),nrow=n)
solve(t(X)%*%X)%*%t(X)
```
```{r}
n <- 3
m <- 3 # m = n
X <- matrix(rnorm(n*m),nrow=n)
solve(t(X)%*%X)%*%t(X)
```
```{r}
n <- 2
m <- 3
X <- matrix(rnorm(n*m),nrow=n)
# solve(t(X)%*%X) %*% t(X) # Error
tX <- t(X)
t(solve(t(tX)%*%tX)%*%t(tX))
t(X)%*%solve(X%*%t(X))
```

Rでは、MASS パッケージにginv() 関数があり、Xの行数・列数に関係なく、ムーアペンローズ疑似逆行列を算出できる。

任意の行列を次のように分解することを特異値分解と言うが

$$
X = U\Sigma V^T
$$

ムーアペンローズ疑似逆行列は
$$
gen.inv(X) = V \Sigma^{+} U^T
$$
で与えられる。ただし$\Sigma$は対角成分に特異値を持つ行列であり、$/Simga^{+}$は特異値の逆数を対角成分に持つ。

実際、Rではこれを利用して、MASSパッケージのginv()関数(g:generalized, inv:inverse)がXの行数・列数を場合分けすることなく、ムーアペンローズ疑似逆行列を返す。

```{r}
library(MASS)
n <- 3
m <- 2 # m < n
X <- matrix(rnorm(n*m),nrow=n)
solve(t(X)%*%X)%*%t(X)
ginv(X)
```
```{r}
n <- 3
m <- 3
X <- matrix(rnorm(n*m),nrow=n)
solve(t(X)%*%X)%*%t(X)
ginv(X)
```
```{r}
n <- 2
m <- 3
X <- matrix(rnorm(n*m),nrow=n)
tX <- t(X)
t(solve(t(tX)%*%tX)%*%t(tX))
ginv(X)
```


# Exercise 1

ムーアペンローズ疑似逆行列の解の幾何的な意味を以下の手順で確認せよ。

## Exercise 1-1

$n=m$

$$
\begin{pmatrix}y_1\\y_2\end{pmatrix} = \begin{pmatrix}x_{11},x_{12}\\x_{21},x_{22}\end{pmatrix}\begin{pmatrix}a_1\\a_2\end{pmatrix}
$$

今、$a_1,a_2$の値を求めるというのは、$(a_1,a_2)$を次元座標として、2次元平面にある2本の直線の交点座標を求めることである。
$$
x_{11} a_1 + x_{12}a_2=y_1\\
x_{21} a_1 + x_{22}a_2=y_2
$$

今、$y_1=3,y_2=1,x_{11}=4,x_{12}=2,x_{21}=-1,x_{22}=3$のとき
の2直線を描き、その交点をMASS::ginv()関数を求め、その点$(\hat{a_1},\hat{a_2})$が交点にあることを、打点することによって示せ。

## Exercise 1-2

$n > m$
$$
\begin{pmatrix}y_1\\y_2\\y_3\end{pmatrix} = \begin{pmatrix}x_{11},x_{12}\\x_{21},x_{22}\\x_{31},x_{32}\end{pmatrix}\begin{pmatrix}a_1\\a_2\\end{pmatrix}
$$

これは3直線の場合である。

$\begin{pmatrix}y_1\\y_2\\y_3\end{pmatrix} = \begin{pmatrix}1\\2\\3\end{pmatrix}$,$\begin{pmatrix}x_{11},x_{12}\\x_{21},x_{22}\\x_{31},x_{32}\end{pmatrix} = \begin{pmatrix}2,3\\1,-2\\-3,-4\end{pmatrix}$とする。


3直線を描図し、ムーアペンローズ疑似逆行列による解$(\hat{a_1},\hat{a_2})$を打点せよ。

## Exercise 1-3

Excerise 1-3 の例で、第1の直線は$(a_1,a_2)$座標について、$y_1 = x_{11}a_1 + x_{12}a_2$で表されているのに対し、解$(\hat{a_1},\hat{a_2})$は、別の直線、$\hat{y_1} = x_{11}\hat{a_1} + x_{12}\hat{a_2}$を表している。

2本の直線を描け。

同様に、第1、第2、第3の直線と、それに対応する、$(\hat{a_1},\hat{a_2})$を通る3直線を描け。

## Exercise 1-4

第1の直線は、$(a_1 = \frac{y_1}{x_{11}},0)$と$(0,\frac{y_1}{x_{12}})$を通る直線である。

それに対応する$(\hat{a_1},\hat{a_2})$を通る直線は$(a_1 = \frac{\hat{y_1}}{x_{11}}=\frac{\hat{a_1}}{x_{11}},0)$と$(0,\frac{\hat{y_1}}{x_{12}}=\frac{\hat{a_2}}{x_{12}})$を通る直線である。
これらは平行である。

点$(\hat{a_1},\hat{a_2})$を通り、この2直線と垂直な直線を引け。その直線の傾きは、$(x_{12},-x_{11})$である。
この垂直な直線が第1の直線と交わる交点$P_1$を求めよ。

$(\hat{a_1},\hat{a_2})$と交点$P_1$との距離を求めよ。

同様に、行列が定める3つの直線のそれぞれについて、行え。

$||y_i,\hat{y_i}||$とこの距離の関係は何か。

## Exercise 1-5

$(a_1,a_2)$平面の任意の点に関して、その点から、直線への垂線の足を求める関数を作成せよ。

その関数を用いて、行列が定める3直線への足、3点が求まる。
$||y_i,\hat{y_i}||$がわかるので、$\sum_{i=1}^3 (y_i-\hat{y_i})^2$も求まるはずである。

この値が$(a_1,a_2)$平面に、どのような高低を作るかを図示し、$(\hat{a_1},\hat{a_2})$の意味を確認せよ。

## Exercise 1-6

$n < m$
$$
\begin{pmatrix}y_1\end{pmatrix} = \begin{pmatrix}x_{11},x_{12} \end{pmatrix}\begin{pmatrix}a_1\\a_2\end{pmatrix}
$$
この場合、直線は1本引ける。
ムーアペンローズ疑似逆行列の解は直線上の1点である。

$||(a_1,a_2)||$が原点から最短距離になっていることを確かめよ。

## Exercise 1-7

$n=m$,$n>m$,$n<m$の3つの場合について、ムーアペンローズ疑似逆行列の解の幾何学的な意味を説明せよ。

## Exercise 1-8

$n<m$の場合に、$y=X\mathbf{a}$の解は、直線上のいずれの点でもよかったが、ムーアペンローズ疑似逆行列は、原点からの距離が最短になるような点$||(\hat{a_1},\hat{a_2})||$を選んだ。

正規化手法と呼ばれる手法では、
$$
\sum_{i=1}^n (y_i-\hat{y_i})^2 + \lambda ||\mathbf{a}||^p
$$
を最小にすることで、残差$\sum_{i=1}^n (y_i-\hat{y_i})^2$を小さくことのみを目指すのではなく、$||\mathbf{a}||^p$という形で、解の取り方に制約を持たせること、その制約の強さをパラメタ$\lambda$でコントロールする。

ムーアペンローズによる解の選び方と、正則化手法での解の選び方との類似性についてコメントせよ。

2016-12-26 統計遺伝学分野で生き抜くための数学〜精選6〜

[][][][]行列でPCA

---
title: "PCA with Matrix 行列でPCA"
author: "ryamada"
date: "2016年12月24日"
output: 
  html_document:
    toc: true
    toc_depth: 6
    number_section: true
---

```{r setup, include=FALSE}
library(knitr)
opts_chunk$set(echo = TRUE)
library(rgl)
knit_hooks$set(rgl = hook_webgl)
```

# 行列の特徴

正方行列の特徴を思い出そう。

正方行列は回転させても変わらない性質があった。

固有値の和(トレース)と固有値の積(行列式)がそれである。

『行列の特徴』のExercises 4-6, 4-7で見たように、分散共分散行列を回転させると、固有値の和と積は変わらないが、固有値の内訳は変わる。
そして、対角化すると、固有値の大小が最も大きくばらつくようにできることを見た。
最も大きくばらつくとは、最大固有値が、色々な回転の中で最大となり、最小固有値が色々な回転の中で最小になるような内訳の取り方を指す。

#Exercise 1
## Exercise 1-1 
『行列の特徴』のExercises 4-6, 4-7を再度実施せよ。

# PCA

PCAは分散共分散行列を回転させて、固有値のばらつきを大きくして、数少ない固有値で、すべての固有値の和の多くの部分を説明しようとするものである。


### Exercise 4-6

xの分散共分散行列の固有値を求めることで、分散共分散行列を対角化するような回転をしたときの、各軸の分散を求めることができる。

その分散を4-4のプロットに重ねてプロットせよ。

```{r,echo=FALSE}
x.cov <- cov(x)
lambdas <- sort(eigen(x.cov)[[1]])
matplot(t(sorted.vars),type="l",xlab="order",ylab="var",col=1)
points(1:d,lambdas,col=2,pch=20)
for(i in 1:(d-1)){
  segments(i,lambdas[i],i+1,lambdas[i+1],col=2)

}
```

### Exercise 4-7

4-4,4-5と同様のプロットを変数の数を増やして実施せよ。

```{r}
d <- 10
n <- 1000
x <- matrix(rnorm(d*n),ncol=d)
x <- apply(x,2,cumsum)

plot(as.data.frame(x[,1:5]),pch=20,cex=0.1)
```

```{r echo=FALSE}
n.iter <- 100
sorted.vars <- matrix(0,n.iter,d)
for(i in 1:n.iter){
  R <- Random.Start(d)
  x.tmp <- t(R %*% t(x))
  cov.tmp <- cov(x.tmp)
  sorted.vars[i,] <- sort(diag(cov.tmp))
}
matplot(t(sorted.vars),type="l",xlab="order",ylab="var")
```
```{r,echo=FALSE}
x.cov <- cov(x)
lambdas <- sort(eigen(x.cov)[[1]])
matplot(t(sorted.vars),type="l",xlab="order",ylab="var",col=1,ylim=range(c(sorted.vars,lambdas)))
points(1:d,lambdas,col=2,pch=20)
for(i in 1:(d-1)){
  segments(i,lambdas[i],i+1,lambdas[i+1],col=2)

}
```