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

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

2017-02-18 統計遺伝学分野で生き抜くための数学〜精選14〜

[][][][]ヤコビ行列と座標変換

---
title: "Calculus7 Jacobian Matrix and Coordinate Transformation ヤコビ行列と座標変換"
author: "ryamada"
date: "2017年2月15日"
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)
```

# 多変数関数の傾き ヤコビ行列 ヘッセ行列 Jacobian and Hessian Matrices

$y = f(x)$という1変数(xの数が1)のスカラー関数(yが1変数)を$x$で微分すると、それは、$y=f(x)$という曲線の傾きを表す。

ヤコビ行列とヘッセ行列はその多変数版。

The derivative of a one-variable scaler function $y = f(x)$ indicates its slope.

Jacobian and Hessian Matrices are multi-variable versions of the slope.

# 2つの用途 Two Ways to use

遺伝統計学分野にとってのヤコビ行列は、傾きを表す行列としての意味と、変数変換・座標変換に関する情報を表す行列としての意味との2つの面で重要である。

There are two important ways to use Jacobian matrix in statistical genetics; one is the indicator of slope and the other is for variable transformation/coordinate transformation.

まず、傾きを表す行列としてのヤコビ行列(とヘッセ行列)を扱い、ついで、変数変換・座標変換でのヤコビ行列を扱う。

Its use for slope will be handled first and its use for coordinate transformation will come later.

# 傾きとしてのヤコビ行列 Jacobian Matrix for Slope

## 正規分布パラメタ推定 Parameter Estimation of Normal Distribution

正規分布のパラメタを最尤推定することにする。

平均$m$、SD $s$の正規分布からのサンプル $X=\{x_1,x_2,...,x_n\}$がある。

Assume the task of maximum likelihood estimation of parameters of normal distribution, mean $m$ and sd $s$ for samples $X=\{x_1,x_2,...,x_n\}$.


尤度関数は The likelihood function is;

$$
L(m,s|X) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi}s}e^{-\frac{(x_i-m)^2}{2 s^2}} . 
$$


対数尤度関数は The logarithm of it is;

$$
LL(m,s|X) = n \times  (-\frac{1}{2}\log{2\pi} - \log{s}) - \sum_{i=1}^n \frac{(x_i-m)^2}{2 s^2}\\
=-\frac{1}{2s^2}(n\times m^2 -2\sum_{i=1}^n x_i \times  m + \sum_{i=1}^n x_i^2) -n \log{s} -\frac{n}{2}\log{2\pi}
$$
この尤度関数と対数尤度関数は、$m,s$という2変数スカラー値関数(スカラー値が尤度)。

These functions are scaler functions with two variables $m$ and $s$.

## 具体例 Example

$X= \{0,1,2\}$という観察をしたとする。

Assume $X= \{0,1,2\}$.

$$
LL(m,s|X) = -\frac{1}{2s^2} (3m^2-6m+5) -3\log{s} -\frac{3}{2}\log{2\pi}
$$
$(m,s)$の最尤推定解 $(\hat{m},\hat{s})$では、$LL$の偏微分が0($\frac{\partial LL(\hat{m},\hat{s})}{\partial m} =0,\frac{\partial LL(\hat{m},\hat{s})}{\partial s} =0$)であるから、以下の連立方程式をとけばよい。

We will get the maximum likelihood estimates $(\hat{m},\hat{s})$, by solving the system of two equations that are partial derivatives of $LL(m,s|X)$.

1つの(対数)尤度関数を各変数で偏微分して連立方程式を作ると、変数の数の方程式が得られる。

この例では、2つの変数、2つの方程式。

When we have one scaler (logarithm) likelihood function and partially differentiate it with $n$ parameters, we will have a system of $n$ equations.

In this case, two parameters and the equations.


$$
\frac{\partial LL(m,s)}{\partial m} = -\frac{1}{2s^2}(6m-6) = 0\\
\frac{\partial LL(m,s)}{\partial s} = \frac{1}{s^3}(3m^2-6m+5) - \frac{3}{s} = 0
$$

第1式より$\hat{m}=1$、それを用いて第2式より$\hat{s}^2 = \frac{2}{3}$を得る。

This can be solved easily as $\hat{m}=1$ and $\hat{s}^2 = \frac{2}{3}$.

### ヤコビ行列・ヘッセ行列 Jacobian and Hessian Matrices

この例のように、
$n$個の変数 $x_1,...,x_n$があって、同数$n$個の関数 $f_1(x_1,...,x_n),f_2(x_1,...,x_n),...,f_n(x_1,...,x_n)$ があるとき、

各関数を各変数で偏微分し、その偏微分関数を$n \times n $行列に並べたもの
$$
\partial F = \begin{pmatrix} \frac{\partial}{\partial x_1}f_1, ...,\frac{\partial}{\partial x_1}f_n\\
...,...,...\\
\frac{\partial}{\partial x_n}f_1, ...,\frac{\partial}{\partial x_n}f_n \end{pmatrix} = \begin{pmatrix}\frac{\partial}{\partial x_i}f_j \end{pmatrix}
$$

を$n$個の変数と$n$個の関数のヤコビ行列と言う。

When we have $n$ functions of $n$ variables, we can make $n\times n$ matrix whose elements are partial derivatives of each equation by each variable, as above. This $n\times n$ matrix whose elements are functions is called Jacobian matrix for the system of equations.


上の例では、1つの$n=2$変数スカラー関数を各変数で偏微分することによって$n$個の関数を得たが、それを$n$個の変数でされに偏微分することができるが、そのようにしてできた、$n\times n$行列を、元のスカラー関数にとってのヘッセ行列と言う。

We made the system of two equations by partially differentiate one scaler function of two variables and then made its Jacobian matrix. This Jacobian matrix is called Hessian matrix of the initial scaler function.

## 接線・接面 Tangent Line and Tangent Plane

### 2x2 = 4 次元空間 Four dimensional space (2x2 = 4)

2つの方程式があった。

Recall the equations.
$$
f_m(m,s) = -\frac{1}{2s^2}(6m-6) = 0\\
f_s(m,s) = \frac{1}{s^3}(3m^2-6m+5) - \frac{3}{s} = 0
$$



$f_m(m,s) = u,f_s(m,s)=v$と書き直すと、

$m,s$が作る2次元平面の上に$(u,v)$という2次元空間上の点が対応づいて広がっているので、この様子を表示するには$(m,s,u,v)$という4次元空間が必要である。

Rewriting as $f_m(m,s) = u,f_s(m,s)=v$, each point $(m,s)$ in a two dimensional space has a value set $(u,v)$. This corresponds to a point $(m,s,u,v)$ which is a point in 4 dimensional space. It is not possible to visualize the 4 dimensional data as they are but we can draw something in 3D space as below.

少し工夫をして描いてみる。
u,vをそれぞれ別々に描けば
```{r,rgl=TRUE}
M <- seq(from=-1,to=2,length=80)
S <- seq(from=0.5,to=1.5,length=80)

ms <- expand.grid(M,S)
m <- ms[,1]
s <- ms[,2]
u <- -1/(2*s^2)*(6*m-6)
v <- 1/(s^3)*(3*m^2-6*m+5)-3/s
plot3d(m,s,u,main="u over m-s plane")
```
```{r rgl=TRUE}
plot3d(m,s,v,main = "v over m-s plane")
```

s=1に固定して、mと$f_m=u,f_s=v$とを描くと

Fix s at 1 and draw (u,v) along m axis,

```{r,rgl=TRUE}
s1 <- 1
us1 <- -1/(2*s1^2)*(6*M-6)
vs1 <- 1/s1^3*(3*M^2-6*M+5)-3/s1
plot3d(M,us1,vs1,type="l",main="(u,v) along m when s=1")
```

m=0.8に固定して、sと$f_m=u,f_s=v$とを描くと

Fix m at 0.8 and draw (u,v) for s axis,
```{r,rgl=TRUE}
m0 <- 0.8
um0 <- -1/(2*S^2)*(6*m0-6)
vm0 <- 1/S^3*(3*m0^2-6*m0+5)-3/S
plot3d(S,um0,vm0,type="l",main="(u,v) along s when m=0.8")
```

### 接面 Tangents


点$(m,s)$において、$(u=f_m(m,s),v=f_s(m,s))$ がどのように傾いているのかを調べる。

Let's evaluate tangents at $(m,s,u=f_m(m,s),v=f_s(m,s))$.


$s$を固定し、$m$方向での$u=f_m$の変化と$v=f_s$の変化を評価すると、それぞれ

Fix s and evaluate change of $u=f_m$ and $v=f_s$ along $m$;

$$
\frac{\partial}{\partial m}f_m(m,s) = -\frac{3}{s^2}\\
\frac{\partial}{\partial m}f_s(m,s) = \frac{6(m-1)}{s^3}
$$

$m$を固定し、$s$方向での$u=f_m$の変化と$v=f_s$の変化を評価すると、それぞれ

Fix m and evaluate change of $u=f_m$ and $v=f_s$ along $s$;
$$
\frac{\partial}{\partial s}f_m(m,s) = \frac{1}{s^3}(6m-6)\\
\frac{\partial}{\partial s}f_s(m,s) = -\frac{3}{s^4}(3m^2-6m+5) + \frac{3}{s^2}
$$
であり、結局、

These are summarized as;

$$
\begin{pmatrix} \frac{\partial}{\partial m}f_m(m,s),\frac{\partial}{\partial s}f_m(m,s)\\
\frac{\partial}{\partial m}f_s(m,s), \frac{\partial}{\partial s}f_s(m,s) \end{pmatrix} \\
=\begin{pmatrix} \frac{\partial ^2 }{\partial m\partial m}f(m,s),\frac{\partial^2}{\partial s\partial m}f(m,s)\\
\frac{\partial^2}{\partial m\partial s}f(m,s), \frac{\partial^2}{\partial s\partial s}f(m,s) \end{pmatrix}\\
= \begin{pmatrix} -\frac{3}{s^2}, \frac{1}{s^3}(6m-6) \\ \frac{6(m-1)}{s^3}, -\frac{3}{s^4}(3m^2-6m+5)+\frac{3}{s^2}\end{pmatrix}
$$

となる。

$f_m,f_s$に関するヤコビ行列であり、対数尤度関数に関するヘッセ行列である。

This is the Jacobian matrix for $f_m,f_s$ and the Hessian matrix for the log-likelihood function.

点$(mx,sx)$における接面を描いてみる。

Let's draw the tangents plane(s) at $(mx,sx)$.

```{r,rgl=TRUE}
mx <- m0
sx <- s1
fm.mxsx <- -1/(2*sx^2)*(6*mx-6)
fs.mxsx <- 1/(sx^3) *(3*mx^2-6*mx+5)-3/sx

dfmdm <- -3/(sx^2)
dfmds <- 1/(sx^3)*(6*mx-6)
dfsdm <- 6*(mx-1)/(sx^3)
dfsds <- -3/(sx^4)*(3*mx^2-6*mx+5)+3/(sx^2)
```
```{r rgl=TRUE}
plot3d(ms[,1],ms[,2],u)
spheres3d(mx,sx,fm.mxsx,radius = 1,col=3)
u.tan <- fm.mxsx + dfmdm * (m-mx) + dfmds * (s-sx)
spheres3d(m,s,u.tan,radius = 0.5,col=2)
```
```{r,rgl=TRUE}
plot3d(m,s,v)
spheres3d(mx,sx,fs.mxsx,radius = 10,col=3)
v.tan <- fs.mxsx + dfsdm * (m-mx) + dfsds * (s-sx)
spheres3d(m,s,v.tan,radius = 5,col=2)
```

s=1に固定して、mと$f_m=u,f_s=v$とを描くと

Fix s at 1 and draw (u,v) along m;
```{r,rgl=TRUE}
plot3d(M,us1,vs1,type="l")
spheres3d(mx,fm.mxsx,fs.mxsx,radius = 1,col=3)
us1.tan <- fm.mxsx + dfmdm * (M-mx) + dfmds * (s1-sx)
vs1.tan <- fs.mxsx + dfsdm * (M-mx) + dfsds * (s1-sx)
spheres3d(M,us1.tan,vs1.tan,radius = 0.1,col=2)
```



m=0に固定して、sと$f_m=u,f_s=v$とを描くと

Fix m at 0 and draw (u,v) along s;
```{r,rgl=TRUE}
plot3d(S,um0,vm0,type="l")
spheres3d(sx,fm.mxsx,fs.mxsx,radius = 0.5,col=3)
um0.tan <- fm.mxsx + dfmdm * (m0-mx) + dfmds * (S-sx)
vm0.tan <- fs.mxsx + dfsdm * (m0-mx) + dfsds * (S-sx)
spheres3d(S,um0.tan,vm0.tan,radius = 0.1,col=2)
```


## ニュートン法 Newton's method

上記の正規分布変数推定の例では、連立方程式は、簡単に解けてしまったが、ニュートン法を用いて、数値計算的に解を求める練習をすることにする。

The problem to estimate parameters of normal distribution was easy and we already know the answer. The followings are its numeric solution based on the Newton method as practice.


1変数関数 $f(x)=0$のニュートン法では、適当なxを初期値と定め、そこの接線を取り、接線が0と交わる値を次のxの値として取り、それを繰り返した。

In case of the Newton method for one parameter scaler function, we set an initial value and drew the tangent line then renewed the value with the coordinate where the tangent line crosses the line $0$ and repeated the procedures.


今回は(m,s)2変数があり、2つの式があり、$(f_m(m,s),f_s(m,s))=(0,0)$となるように$(m,s)$の値の組を更新したい。

In our case we have two parameters to be estimated and we have two functions to handle $(f_m(m,s),f_s(m,s))$. We want to repeat procedures to renew $(m,s)$.


正規分布のパラメタ推定は、式がややこしいので、簡易な例をとることにする。

Here we will take a simpler function rather than the complicated function from normal distribution example.

最小化したい2変数スカラー関数を考える。

The following is the function.

$$
f(x,y) = (x+y)^4 + (x-y)^4
$$
```{r}
fr <- function(w) {
    x <- w[1]
    y <- w[2]
    (x+y)^4 + (x-y)^4
}
t <- seq(from=-2,to=2,length=50)
xy <- as.matrix(expand.grid(t,t))
x <- xy[,1]
y <- xy[,2]
z <- apply(xy,1,fr)
plot3d(x,y,z)
```



答えはわかっていて、$x=y=0$。

The answer is known as $x=y=0$.

それを偏微分を用いて計算機的に求めてみる。

Let's solve it numerically with partial derivatives.


$$
\frac{\partial f}{\partial x} = f_x = 4(x+y)^3 + 4(x-y)^3\\
\frac{\partial f}{\partial y} = f_y = 4(x+y)^3 - 4(x-y)^3
$$

最小値をとる$(\hat{x},\hat{y})$では、$f_x(\hat{x},\hat{y})= 0, f_y(\hat{x},\hat{y})=0$。

ニュートン法を実施するために
$f_x,f_y$をさらに偏微分する。

Partially differentiate $f_x,f_y$ further with x and y.


$$
\frac{\partial f_x}{\partial x} = \frac{\partial ^2 f}{\partial x\partial x} = f_{xx} = 12(x+y)^2 + 12 (x-y)^2\\
\frac{\partial f_x}{\partial y} = \frac{\partial ^2 f}{\partial x\partial y} = f_{xy} = 12(x+y)^2 - 12 (x-y)^2\\
\frac{\partial f_y}{\partial x} = \frac{\partial ^2 f}{\partial y\partial x} = f_{yx} = 12(x+y)^2 - 12 (x-y)^2\\
\frac{\partial f_y}{\partial y} = \frac{\partial ^2 f}{\partial y\partial y} = f_{yy} = 12(x+y)^2 + 12 (x-y)^2
$$

$(x_0,y_0,f_x(x_0,y_0),f_y(x_0,y_0))$における接面と$f_x=0,f_y=0$との交点座標を$(x_1,y_1,0,0)$とする。

The tangent plane at $(x_0,y_0,f_x(x_0,y_0),f_y(x_0,y_0))$ crosses with the plane $f_x=0,f_y=0$; the crossing point is now called as $(x_1,y_1,0,0)$ and $(x_1,y_1)$ should replace $(x_0,y_0)$.


接面が$f_1=f_2=0$をよぎる$(m,s)$の座標を取り出すことは、

$(x_1,y_1)$ can be obtained by soving the equations below;

$$
f_x(x_1,y_1) = f_x(x_0,y_0) + \frac{\partial f_x}{\partial x} (x_1-x_0) + \frac{\partial f_x}{\partial y} (y_1-y_0) = 0\\
f_y(x_1,y_1) = f_y(x_0,y_0) + \frac{\partial f_y}{\partial x} (x_1-x_0) + \frac{\partial f_y}{\partial y} (y_1-y_0) = 0\\
$$

を解くことであるから、

行列
$$
\partial F = \begin{pmatrix} \frac{\partial}{\partial x}f_x(x,y),\frac{\partial}{\partial y}f_x(x,y)\\
\frac{\partial}{\partial y}f_y(x,y), \frac{\partial}{\partial y}f_y(x,y) \end{pmatrix} = \begin{pmatrix} f_{xx},f_{xy}\\ f_{yx},f_{yy}\end{pmatrix} = 12 \begin{pmatrix} (x+y)^2+(x-y)^2,(x+y)^2-(x-y)^2\\(x+y)^2-(x-y)^2,(x+y)^2+(x-y)^2 \end{pmatrix}
$$



を用いると、

Using this matrix, the Hessian matrix, the system is expressed as below;

$$
\begin{pmatrix} f_x(x_0,y_0) \\ f_y(x_0,y_0)\end{pmatrix} + \partial F \begin{pmatrix} x_1- x_0 \\ y_1- y_0 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}
$$

と表せる。

この行列 $\partial F$は式$f(x,y) = (x+y)^4 + (x-y)^4$のヘッセ行列である。



$$
\begin{pmatrix} x\\y \end{pmatrix} = w
$$
とすれば

Using this notation, the equation is expressed as below.

$$
F(w_0) + \partial F (w_1-w_0) = 0
$$

となり、結局
$$
\partial F w_1 = \partial F w_0 - F(w_0)\\
w_1 = w_0 - (\partial F)^{-1} F(w_0)
$$
によって求まる。

たとえば、

The following is an example.

$$
(x_0,y_0) = (3,4)
$$
のとき
$$
\partial F(3,4) =  12 \begin{pmatrix} (x+y)^2+(x-y)^2,(x+y)^2-(x-y)^2\\(x+y)^2-(x-y)^2,(x+y)^2+(x-y)^2 \end{pmatrix} = 12 \begin{pmatrix} 50, 48\\48,50\end{pmatrix}
$$
であるから、
$$
(\partial F(3,4))^{-1} = (12 \begin{pmatrix} 50, 48\\48,50\end{pmatrix})^{-1} .
$$

実際に計算すれば

The inverse matrix can be calculated;

```{r}
fx.F <- function(x,y){
  tmp1 <- (x+y)^2
  tmp2 <- (x-y)^2
  12 * matrix(c(tmp1+tmp2,tmp1-tmp2,tmp1-tmp2,tmp1+tmp2),2,2)
}
x0 <- 1
y0 <- 2
F.inv <- solve(fx.F(x0,y0))
F.inv
```


となり、

```{r}
w0 <- matrix(c(x0,y0),nrow=2)

fxfy <- function(w){
  c(4*(w[1]+w[2])^3+4*(w[1]-w[2])^3,4*(w[1]+w[2])^3-4*(w[1]-w[2])^3)
}

w1 = w0 - F.inv %*% matrix(fxfy(w0),nrow=2)
w1
```


この$w_1=(x_1,y_1)$を使って$(x_1,y_1,f_x(x_1,y_1),f_y(x_1,y_1))$の接面と$f_x=0,f_y=0$との交点座標を$(x_2,y_2,0,0)$とすれば、

The $(x,y)$ can be renewed further.

```{r}
w2 <- w1 - solve(fx.F(w1[1],w1[2])) %*% matrix(fxfy(w1),nrow=2)
w2
```


これを繰り返すことにする。

Repeat the procedures as below.

```{r}
x0 <- 3
y0 <- 4
w <- matrix(c(x0,y0),nrow=2)
xy <- matrix(c(x0,y0),nrow=1)
n.iter <- 10
for(i in 1:n.iter){
  new.w <- w - solve(fx.F(w[1],w[2])) %*% matrix(fxfy(w),nrow=2)
  w <- new.w
  xy <- rbind(xy,c(w))
}
plot(xy,type="b",pch=20)
matplot(xy,type="l")
```

## Exercise 1

### Exercise 1-1

正規分布の変数$m,s$にニュートン法を適用してみる。以下の説明、コードを理解せよ。

The followings are Newton's method applied to the problem of normal distribution. Read and understand the steps and the codes.


正規分布の変数推定対数尤度関数から、2つの方程式$f_m(m,s)=0$, $f_s(m,s) = 0$が得られる。

Two equations were obtained from Log-likelihood function.

このヤコビ行列を作る。

Make the Jacobian matrix.

ヤコビ行列を利用して、
接面が$u=f_1(m,s)=0$, $v=f_2(m,s)=0$と交わる$(m,s)$を採用する。

Renew the initial (m,s) using the Jacobian matrix.



```{r rgl=TRUE}
mx <- 0.6
sx <- 1
plot3d(m,s,u,col=gray(0.3),size=0.01)
spheres3d(mx,sx,fm.mxsx,radius = 1,col=3)
u.tan <- fm.mxsx + dfmdm * (m-mx) + dfmds * (s-sx)
spheres3d(m,s,u.tan,radius = 0.5,col=2)
spheres3d(m,s,rep(0,length(m)),radius = 0.5,col=4)
```


接面が$u=0$平面と交わって直線が現れる。

The tangent plane corsses with the plane $u=0$, showing a straight line.


```{r,rgl=TRUE}
plot3d(m,s,v,col=gray(0.3),size=0.01)
spheres3d(mx,sx,fs.mxsx,radius = 10,col=3)
v.tan <- fs.mxsx + dfsdm * (m-mx) + dfsds * (s-sx)
spheres3d(m,s,v.tan,radius = 1,col=2)
spheres3d(m,s,rep(0,length(m)),radius = 1,col=4)

```


接面が$v=0$平面と交わって直線が現れる。

The tangent plane corsses with the plane $v=0$, showing a straight line.

この2つの直線の交点が、ニュートン法が示す、次の$(m1,s1)$の位置。

$(m1,s1)$ is the crossing point of these two lines.

```{r,rgl=TRUE}
plot3d(mx,sx,fm.mxsx,radius = 1,col=3,zlim = c(-10,10))
spheres3d(mx,sx,fs.mxsx,radius = 1,col=5,zlim = c(-10,10))
spheres3d(m,s,u.tan,radius = 0.2,col=3,zlim = c(-10,10))
spheres3d(m,s,v.tan,radius = 0.2,col=5,zlim = c(-10,10))
spheres3d(m,s,rep(0,length(m)),radius = 0.2,col=4,zlim = c(-10,10))

```

```{r}
matrix(c(dfmdm,dfmds,dfsdm,dfsds),byrow=TRUE,nrow=2)
```

```{r}
new.ms <- matrix(c(mx,sx),nrow=2) - solve(matrix(c(dfmdm,dfmds,dfsdm,dfsds),byrow=TRUE,nrow=2)) %*% matrix(c(fm.mxsx,fs.mxsx),nrow=2)
new.ms
```



```{r,rgl=TRUE}
plot3d(mx,sx,fm.mxsx,radius = 0.5,col=3)
spheres3d(mx,sx,fs.mxsx,radius = 0.5,col=5)
spheres3d(m,s,u.tan,radius = 0.2,col=3)
spheres3d(m,s,v.tan,radius = 0.2,col=5)
spheres3d(m,s,rep(0,length(m)),radius = 0.2,col=1)
spheres3d(new.ms[1],new.ms[2,],0,radius = 0.5,col=6)
```

$\partial F$を計算する関数を使って
```{r}
partialF <- function(m,s){
  matrix(c(-3/(s^2),1/(s^3)*(6*m-6),6*(m-1)/(s^3),-3/(s^4)*(3*m^2-6*m+5) + 3/(s^2)),byrow=TRUE,2,2)
}
m0 <- mx
s0 <- sx
n.iter <- 10
ms <- matrix(c(m0,s0),nrow=1)
solve(partialF(ms[1,1],ms[1,2]))

for(i in 1:n.iter){
  pF <- partialF(ms[i,1],ms[i,2])
  tmpf1.mxsx <- -1/(2*ms[i,2]^2)*(6*ms[i,1]-6)
  tmpf2.mxsx <- 1/(ms[i,2]^3) *(3*ms[i,1]^2-6*ms[i,1]+5)-3/ms[i,2]
  new.ms <- matrix(ms[i,],nrow=2) - solve(pF) %*% matrix(c(tmpf1.mxsx,tmpf2.mxsx),nrow=2)
  ms <- rbind(ms,c(new.ms))
}

plot(ms,type="b")
points(1,sqrt(2/3),pch=20,col=2)
matplot(ms,type="l")
abline(h=c(1,sqrt(2/3)),col=3)
```

# 座標変換とヤコビ行列 Coordinate Transformations and Jacobian Matrix

$$
u = \frac{\partial}{\partial x} f(x,y) = f_{x}(x,y)\\
v = \frac{\partial}{\partial x} f(x,y) = f_{y}(x,y)
$$
この関係を

These equations can be considered as coordinate transformation.

$$
u = f(x,y)\\
v = g(x,y)
$$
と捉えなおす。

逆に The opposite transformation;
$$
x = h(u,v)\\
y = k(u,v)
$$
も成り立つとする。

先ほどまで$(x,y)$と言う座標で呼ばれてた点が、$(u,v)$という座標で呼ばれるようになった。言い換えれば、そのような座標を引きなおした、とも読み取れる。

The point who had $(x,y)$ coordinate, has now different coordinate $(u,v)$. On the same space, new coordinate-system lattice appeared, in other words.

$$
u = 2x+\frac{1}{2}y\\
v = x+y
$$

逆に(u,v)という座標でよばれていた点を(x,y)という座標で呼ぶように変えるのが、以下の処理である。

The opposite relation is;

$$
x = \frac{2u-v}{3}\\
y = \frac{-2u + 4v}{3}
$$

## 格子 Lattice

上記の関係だと、新たな座標系による格子は次のようになる。

The lattice with new coordinate system is drawn as below.

左のパネルでは$x$の値が等しい点を結んだ線と、$y$の値が等しい点を結んだ線とが現れている。
右のパネルでは$u$の値が等しい点を結んだ線と、$v$の値が等しい点を結んだ線が現れている。

In the left panel, lines indicate sets of points with same x values or same y values.

In the right panel, lines indicate sets of points with same u values or same v values.

元の座標で点を打ち、元の座標で格子を描くと、マス目状の格子が描かれ、元の座標で点を打ち、新たな座標で格子を描くと、ゆがんだ格子が描かれる。

When points are placed in the original (x,y) coordinates and lattice of the original coordinate, then the lattice appears as perpendicular lattice.

When lattice is drawn with (u,v) system, then the lattice lines are curved.


```{r}
n <- 10^4
XY <- matrix(runif(n*2)*5,ncol=2)
x <- XY[,1]
y <- XY[,2]

col.ori <- abs(floor(x-min(x)+1)*floor(y-min(y)+1))
par(mfcol=c(1,2))
plot(x,y,col=col.ori)
#u <-(x^2-y)/10
#v <- exp((x+y)/10)
u <- 2*x+1/2*y
v <- x+y
col <- abs(floor(u-min(u)+1) * floor(v-min(v)+1))
plot(x,y,col=col,pch=20)
par(mfcol=c(1,1))
```

逆に、新たな座標で点を打ち、それに元の座標で格子を描くと、ゆがみ方が逆になったパターンが現れる。

The opposite;

```{r}
n <- 10^4
UV <- matrix(runif(n*2)*5,ncol=2)
u <- UV[,1]
v <- UV[,2]

col.ori <- abs(floor(u-min(u)+1)*floor(v-min(v)+1))
par(mfcol=c(1,2))
plot(u,v,col=col.ori)
#u <-(x^2-y)/10
#v <- exp((x+y)/10)
x <- (2*u-v)/3
y <- (-2*u+4*v)/3
col <- abs(floor(x-min(x)+1) * floor(y-min(y)+1))
plot(u,v,col=col,pch=20)
par(mfcol=c(1,1))
```

## 積分 Integration

$0 \le x \le 1, 0 \le y \le 1$という正方領域に
$f(x,y) = x^2+y$という関数があって

Assume a function $f(x,y) = x^2+y$ in the area $0 \le x \le 1, 0 \le y \le 1$, that is a square and we want to know the integral below.

$$
\int_{x=0}^{x=1} \int_{y=0}^{y=1} f(x,y) dxdy
$$
という積分をしたいとする。

正方領域を十分に小さな正方形に分割して、足し合わせることで近似してみる。
小さな正方形の面積で重みづけをした和になる。


Divide the square into small squares and sum the value of the points multiplied by their area, then it is good approximation of the integral.

```{r}
n <- 10
tmp <- seq(from=0,to=1,length=n+1)
tmp <- tmp[-(n+1)]
xy <- as.matrix(expand.grid(tmp,tmp))
x <- xy[,1]
y <- xy[,2]
s <- (tmp[2]-tmp[1])*(tmp[2]-tmp[1])
fxy <- x^2+y
sum(fxy * s)
```

分割を細かくすると、真の値に近づくようだ。

With the division finer, the approximation gets closer to the true value.
```{r}
ns <- 2^(1:10)
ret <- rep(0,length(ns))
for(i in 1:length(ns)){
  n <- ns[i]
  tmp <- seq(from=0,to=1,length=n+1)
  tmp <- tmp[-(n+1)]
  xy <- as.matrix(expand.grid(tmp,tmp))
  x <- xy[,1]
  y <- xy[,2]
  deltax <- deltay <- tmp[2]-tmp[1]
  s <- deltax * deltay
  fxy <- x^2+y
  ret[i] <- sum(fxy * s)
}
plot(ns,ret)
```

ここで、小さい面積を$\delta x \times \delta y$として計算した。

We used $\delta x \times \delta y$ as the area of small square.


これは、横軸方向の小さいベクトル$(\delta x,0)$と縦軸方向の小さいベクトル$(0,\delta y)$との作る小さな正方形の面積が$\delta x \times \delta y$と計算できることを利用している。

This squre is consisted of two mutually perpendicular small vectors, $(\delta x,0)$ and $(0,\delta y)$.

では、同様に、$(u,v)$を使うとどうなるだろうか。

Now we want to do the similar for (u,v).

$0 \le u \le 1, 0 \le v \le 1$の領域について、(u,v)の細かい格子を作って
$f(x,y) =f((2*u-v)/3, (-2*u+4*v)/3) = ((2*u-v)/3)^2 +  (-2*u+4*v)/3$ の値を
足し合わせればよい。

Make fine lattice and sum up the values $f(x,y) =f((2*u-v)/3, (-2*u+4*v)/3) = ((2*u-v)/3)^2 +  (-2*u+4*v)/3$ multiplied by their area.

このときの小さい格子の作る小さい領域は
$\delta u = 2 \times \delta x + \frac{1}{2} \delta y$ と言う小さいベクトルと
$\delta v = \delta x - \delta y$と言う小さいベクトルが作る小領域である。

The small area is approxiameted as paralleogram between $\delta u = 2 \times \delta x + \frac{1}{2} \delta y$ and $\delta v = \delta x - \delta y$.


この面積は$|\delta| u \times |\delta v|$というように小ベクトルの長さの積にはならない。
この小面積は平行四辺形なので、平行四辺形の面積の求め方を確認することにする。

The area of parallelogram is not $|\delta| u \times |\delta v|$.

We should know how to calculate the area of parallelogram.



## 平行四辺形の面積 Area lf Parallelogram

ベクトル$A=(a,b)$とベクトル$X=(x,y)$とが作る平行四辺形の面積は、2つのベクトルの成す角を$\theta$とすれば

Assume two vecors $A=(a,b)$ and $X=(x,y)$, the angle between which is $\theta$, that make a parallelogram.

Its area is;
$$
S = |A| |X| \sin{\theta} \\
= |A||X|\sqrt{1-\cos^2{\theta}} \\
= \sqrt{|A|^2|X|^2-(A \cdot X)^2}
$$
ただし、$A\cdot X$は内積。
$(A \cdot X)$ is inner product.


これを成分表示すれば Express the area with elements of vectors;

$$
S = \sqrt{(a^2+b^2)\times (x^2+y^2) - (ax+by)^2}\\
= \sqrt{a^2y^2+b^2x^2 - 2abxy}
= |ay-bx|
$$
この$|ay-bx|$は、行列
$$
\begin{pmatrix}a,x\\b,y\end{pmatrix}
$$
の行列式の絶対値。

$|ay-bx|$ is the absolute value of determinant of the matrix above.

## 変数変換の微小面積と行列式 Small Area after Coordinate Transformation and Determinant of Matrix

変数変換によって求めたい微小面積は
$\delta u = 2 \times \delta x + \frac{1}{2} \delta y$ と言う小さいベクトルと
$\delta v = \delta x + \delta y$とが作る平行四辺形であるから、

The small area for (u,v) coordinate is parallelogram by $\delta u = 2 \times \delta x + \frac{1}{2} \delta y$ and
$\delta v = \delta x + \delta y$.

それは、
$$
\begin{pmatrix} 2 \delta x, \delta x \\ \frac{1}{2} \delta y, \delta y\end{pmatrix} 
$$
という行列の行列式の絶対値。

It is the absolute value of determinant of the matrix above.

それは
$$
2\times \delta x \times \delta y - \delta x \times \frac{1}{2} \delta y = (2\times 1 - 1\times \frac{1}{2}) \times \delta x \delta y = \frac{3}{2} \delta x \delta y.
$$
ここで$2\times 1  - 1 \times \frac{1}{2}$というのは

$$
\begin{pmatrix} 2,1\\\frac{1}{2},1 \end{pmatrix} = \begin{pmatrix} \frac{\partial x}{\partial u} , \frac{\partial x}{\partial v}\\ \frac{\partial y}{\partial x},\frac{\partial y}{\partial v} \end{pmatrix}
$$
の行列式。

これは、ヤコビ行列の行列式であることを意味している。

The above formulas have shown, the matrix for the area calculation is Jacobian matrix of coordinate transformation equations.


結局、

$$
dS = |det(Jacobi)| dudv = dxdy
$$

であるから、

$0 \le u \le 1, 0 \le v \le 1$の領域について、
$f(x,y) =f((2u-v)/3, (-2u+4v)/3) = ((2u-v)/3)^2 +  (-2u+4v)/3$を積分するには
以下を計算することとなる。

The integral of the problem is expressed as;

$$
\int_{u=0}^{u=1} \int_{v=0}^{v=1} ((2u-v)/3)^2 +  (-2u+4v)/3 \times |det \begin{pmatrix} \frac{\partial x}{\partial u} , \frac{\partial x}{\partial v}\\ \frac{\partial y}{\partial x},\frac{\partial y}{\partial v} \end{pmatrix})| dudv \\
= \int_{u=0}^{u=1} \int_{v=0}^{v=1} ((2u-v)/3)^2 +  (-2u+4v)/3 \times |det\begin{pmatrix} 2,1\\\frac{1}{2},-1 \end{pmatrix} | dudv
$$

ここで、
$\begin{pmatrix} 2,1\\\frac{1}{2},-1 \end{pmatrix}$は、u,vの関数ではないので、$dS = C dudv = dxdy$のように、微小面積は座標によらず定数倍($C=\frac{3}{2}$倍)であることがわかるが、それは、図を見ると、平行四辺形の格子がどこも同じ大きさであることに対応している。

In this case, the determinant of Jacobian matrix is constant regardless of (u,v), that were shown in the figure where all parallelograms of lattice were identical.

微小面積の変換が定数倍ではない場合、ヤコビ行列の行列式がu,vの関数である例を練習問題として取り上げる。

The following exercise is the case of non-constant Jacobian determinant.



## Exercise 2
### Exercise 2-1
平面の極座標変換
$$
x = r \cos{\theta}\\
y = r \sin{\theta}
$$
のヤコビ行列を求めよ。

This indicates polar coordinate transformation in 2D. 

Answer its Jacobian matrix.


### Exercise 2-2 
その行列式が $r$であることを示せ。

Show its determinant is $r$.

### Exercise 2-3
半径1の球の体積の半分を求めることを考える。

We want to know the volume of unit sphere.

$x=r\cos{\theta},y=r\sin{\theta}$における、単位球面のz座標は$\pm \sqrt{1-r^2}$であるから、

z-coordinate for $x=r\cos{\theta},y=r\sin{\theta}$ is $\pm \sqrt{1-r^2}$.


半径rを0から1まで、角$\theta$を0から$2\pi$まで$\sqrt{1-r^2}$について積分すれば、単位球の半分の値となる。

Integral $\sqrt{1-r^2}$ from 0 to 1 for r and from 0 to $2\pi$ for $\theta$ is the half of the volume.

微小面積$dS= |det(Jacobi)| drd\theta = r drd\theta$であるから、半球の体積は

The infinitesimal area is $dS= |det(Jacobi)| drd\theta = r drd\theta$, then;

$$
\int_{\theta=0}^{\theta=2\pi} \int_{r=0}^{r=1} \sqrt{1-r^2} (r dr d\theta)
$$
これを用いて、単位球の体積が$\frac{4}{3}\pi$であることを示せ。

Show that the volume of unit sphere is $\frac{4}{3}\pi$, based on the information above.

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

[][][][][]テイラー展開とボンフェロニ補正、積率母関数、ニュートン法

---
title: "Calculus6 Taylor expansion テイラー展開"
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)
```

# 式 Formula

関数が多項式になっている。

A function is expressed as a polynomial function.

$$
\sum_{n=0}^\infty \frac{f^{(n)}(a)}{n!} (x-a)^n
$$

# Examples

$$
f(x) = \sin{x}
$$

```{r}
x <- seq(from=-10,to=10,length=1000)
fx <- sin(x)
plot(x,fx,type="l")
```

$$
f^{(0)}(x) = \sin{x}\\
f^{(1)}(x) = \cos{x}\\
f^{(2)}(x) = - \sin{x}\\
f^{(3)}(x) = - \cos{x}\\
f^{(4)}(x) = \sin{x} = f^{(0)}(x)\\
$$

$$
f^{(4k)}(x) = \sin{x}\\
f^{(4k+1)}(x) = \cos{x}\\
f^{(4k+2)}(x) = - \sin{x}\\
f^{(4k+3)}(x) = - \cos{x}\\
$$
$a=0$のとき、When $a=0$

$$
f^{(4k)}(0) = \sin{0} = 0\\
f^{(4k+1)}(0) = \cos{0} = 1\\
f^{(4k+2)}(0) = - \sin{0} = 0\\
f^{(4k+3)}(0) = - \cos{0} = -1\\
$$
```{r}
fx.0 <- 0/factorial(0) * x^0
plot(x,fx,type="l")
points(x,fx.0,type="l",col=2)
```

```{r}
fx.1 <- fx.0 + 1/factorial(1) * x^1
plot(x,fx,type="l")
points(x,fx.1,type="l",col=2)
```
```{r}
fx.2 <- fx.1 + 0/factorial(2) * x^2
plot(x,fx,type="l")
points(x,fx.2,type="l",col=2)
```
```{r}
fx.3 <- fx.2 + (-1)/factorial(3) * x^3
plot(x,fx,type="l")
points(x,fx.3,type="l",col=2)
```
```{r}
fx.4 <- fx.3 + 0/factorial(4) * x^4
plot(x,fx,type="l")
points(x,fx.4,type="l",col=2)
```
```{r}
fx.5 <- fx.4 + 0/factorial(5) * x^5
plot(x,fx,type="l")
points(x,fx.5,type="l",col=2)
```
```{r}
fx.6 <- fx.5
fx.7 <- fx.6 + 1/factorial(7) * x^7
plot(x,fx,type="l")
points(x,fx.7,type="l",col=2)
```

```{r}
fx.8 <- fx.7
fx.9 <- fx.8 + (-1)/factorial(9) * x^9
plot(x,fx,type="l")
points(x,fx.9,type="l",col=2)
```

# Exercise 1

## Exercise 1-1

Do the same for $a=1$ of $\sin{x}$.

## Exercise 1-2

Do the same for $f(x)= e^x$.

## Exercise 1-3

Do the same of $f(x) = (x-1)(x-2)(x-3)(x-4)$

# Multiple testing correction

$N$個の独立な検定を行い$N$個のp-値を得るとする。
帰無仮説が成り立つとき、すべてのp-値が a以上である確率は

Assume $N$ independent tests. When null hypothesis is true for all, probability that all p-values are equal or more than a is;

$$
f(a) = (1-a)^N
$$
この式の$a=0$周辺の$f(a)^{(1)}$までのテイラー展開は

Taylor expanstion around $a=0$ upto $f(a)^{(1)}$ is;

$$
f(a) \sim \frac{1}{0!}a^0 + (-1)\frac{N}{1!}a^1 = 1- Na
$$

この近似を用いたのがボンフェロニ補正。

Bonferroni's correction is based on this approximation.

# Exercise 2

## Exercise 2-1

以下の関数を$a=0$周辺で$f(a)^k; k=0,1,2,...$までテイラー展開したときの$k$と近似値との関係をプロットし、ボンフェロニ補正が「保守的」であることを説明せよ。

Expand $f(a)$ around $a=0$ upto $f(a)^k; k=0,1,2,...$ and plot the relation between k and the approximated values, then describe that Bonferroni's correction is conservative using the plot.
$$
f(a) = (1-a)^N
$$

# 積率母関数 Moment generating function

## Exercise 3

### Exercise 3-1

このサイト https://www.probabilitycourse.com/chapter6/6_1_3_moment_functions.php を読み、積率母関数とテイラー展開との関係を説明せよ。

Read the site https://www.probabilitycourse.com/chapter6/6_1_3_moment_functions.php .

Describe relation between moment generating function and Taylor expansion.

# ニュートン法 Newton's method

$$
f(x) = \frac{1}{2}e^{x^3} + 2x - 2
$$

```{r}
x <- seq(from=-1.5,to=1.5,length=100)
my.fx <- function(x){
  exp(x^3)/2 + 2*x - 2
}
fx <- my.fx(x)
plot(x,fx,type="l")
abline(h=0,col=2)
```

$f(x)=0$の解を近似的に求める。

Solve $f(x)=0$ numerically.

## 単調性の確認 Monotonicity

# Exercise 4 
## Exercise 4-1 
$f(x)$が単調増加関数であることを示せ。
Show that $f(x)$ is a monotonically increasing function.

## ニュートン法 Newton's method

適当な値$x_0$からスタートし、$\frac{d}{dx}f(x)$を使って、$f(x)=0$の解に次第に近づく。

Take an arbitrary value$x_0$ and get closer to the solution of $f(x)=0$ step-by-step.

```{r}
x0 <- 1

my.dfxdx <- function(x){
  3/2 * x^2 * exp(x^3) + 2
}

this.x <- x0
this.y <- my.fx(this.x)

this.dfxdx <- my.dfxdx(this.x)
this.intercept <- this.y -this.x*this.dfxdx

new.x <- this.x - this.y/this.dfxdx

plot(x,fx,type="l")
points(this.x,this.y,pch=20,col=3)
points(new.x,0,pch=20,col=4)
abline(h=0,col=2)
abline(this.intercept,this.dfxdx,col=3)

x <- seq(from = new.x -0.1,to=this.x+0.1,length=100)
fx <- my.fx(x)
plot(x,fx,type="l")
points(this.x,this.y,pch=20,col=3)
points(new.x,0,pch=20,col=4)
abline(h=0,col=2)
abline(this.intercept,this.dfxdx,col=3)
```

接線とx軸の交点は$f(x)=0$の解に近づいている。

The crossing point with the $y=0$ line of the tangent line at $(x_0,f(x_0))$ gets closer to the crossing point of $y=f(x)$ with the $y=0$ line is closer that $x_0$. 

これを繰り返す。Repeat this procedure.
```{r}
this.x <- new.x
this.y <- my.fx(this.x)

this.dfxdx <- my.dfxdx(this.x)
this.intercept <- this.y -this.x*this.dfxdx

new.x <- this.x - this.y/this.dfxdx

plot(x,fx,type="l")
points(this.x,this.y,pch=20,col=3)
points(new.x,0,pch=20,col=4)
abline(h=0,col=2)
abline(this.intercept,this.dfxdx,col=3)

x <- seq(from = new.x -0.1,to=this.x+0.1,length=100)
fx <- my.fx(x)
plot(x,fx,type="l")
points(this.x,this.y,pch=20,col=3)
points(new.x,0,pch=20,col=4)
abline(h=0,col=2)
abline(this.intercept,this.dfxdx,col=3)
```

```{r}
this.x <- new.x
this.y <- my.fx(this.x)

this.dfxdx <- my.dfxdx(this.x)
this.intercept <- this.y -this.x*this.dfxdx

new.x <- this.x - this.y/this.dfxdx

plot(x,fx,type="l")
points(this.x,this.y,pch=20,col=3)
points(new.x,0,pch=20,col=4)
abline(h=0,col=2)
abline(this.intercept,this.dfxdx,col=3)

x <- seq(from = new.x -0.1,to=this.x+0.1,length=100)
fx <- my.fx(x)
plot(x,fx,type="l")
points(this.x,this.y,pch=20,col=3)
points(new.x,0,pch=20,col=4)
abline(h=0,col=2)
abline(this.intercept,this.dfxdx,col=3)
```

```
}

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

[][][][]曲線・曲面

微分積分5曲線と曲面

微分積分5曲線と曲面

f:id:ryamada22:20170128115550p:image

---
title: "Calculus5 曲線・曲面 Curved line and surface"
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)
```

# 曲線の長さ Length of curve

$$
y = f(x) = x + \frac{1}{2}x^2
$$

```{r}
x <- seq(from=0,to=3,length=100)
y <- x+ 1/2*x^2
plot(x,y,type="l")
```



$$
\frac{d}{dx}f(x) = 1 + x
$$
微分すると、単位x増分当たりのyの増加量がわかる。それらは、微小直角三角形。

Differentiation gives increase of y along with unit increase in x. They are two edges of a right triangle.

三角形 (x,y) ,(x+dx,y),(x+dx,y+dy)は直角三角形の3頂点。
曲線の長さはこの三角形の斜辺の積分

(x,y) ,(x+dx,y),(x+dx,y+dy) are three vertices of a right triangle.
The length of curve is the integral of the obliques.

$$
\int_{0}^3 \sqrt{dx^2+dy^2} = \int_0^3 \sqrt{(\frac{dx}{dx})^2+(\frac{dy}{dx})^2}dx = \sqrt{1+(\frac{dy}{dx})^2}dx \\
$$

積分できるのか?
Can we integrate this?

$$
\int_0^3 \sqrt{2+2x+x^2}dx\\
=[\sqrt{2+2x+x^2}(x+1) + \sinh^{-1}(x+1)]^3_0
$$

Mathematicaを使えば…。

Check Mathematica http://www.wolframalpha.com/examples/Math.html for "Calculus AND ANALYSIS"
and http://d.hatena.ne.jp/ryamada22/20170128 

積分できるかできないかわからなくても近似計算はできる。

Even when the integral is not easy shape, computational calculation is possible.

```{r}
integrate(function(x){(2+2*x+x^2)^(1/2)},0,3)
```

# Exercise 1

## Exercise 1-1

$ y = f(x) = e^{x^2}$ の $x \in [-1,1]$の長さを計算せよ。

Calculate the length of $ y = f(x) = e^{x^2}$ for $x \in [-1,1]$.

# 曲線のパラメタ表示 Parameteric expression of curve.

$$
x = \cos{t}\\
y = \sin{t}\\
t \in [0,2\pi)
$$

```{r}
t <- seq(from=0,to=1,length=100)*2*pi
t <- t[-length(t)]
x <- cos(t)
y <- sin(t)
plot(x,y,type="l")
matplot(cbind(x,y),type="l")
```

$$
L(\tau) = \int_0^\tau \sqrt{(\frac{dx}{dt})^2 + (\frac{dy}{dt})^2} dt
$$

$$
\frac{dx}{dt} = -\sin{t}\\
\frac{dy}{dt} = \cos{t}\\
L(\tau) = \int_0^\tau \sqrt{(\sin{t})^2 + (\cos{t})^2} dt\\
=\int_0^\tau 1 dt\\
=[t]^\tau_0 = \tau
$$

## 同じ軌道、スピードが違う Same trajectory but different speed

$$
x = \cos{t^2}\\
y = \sin{t^2}\\
t \in [0,\sqrt{2\pi})
$$

```{r}
t <- seq(from=0,to=1,length=100)*sqrt(2*pi)
t <- t[-length(t)]
x <- cos(t^2)
y <- sin(t^2)
plot(x,y,type="l")
matplot(cbind(x,y),type="l")
```

$$
L(\tau) = \int_0^\tau \sqrt{(\frac{dx}{dt})^2+(\frac{dy}{dt})^2}dt\\
= \int_0^\tau \sqrt{(-2t\sin(t^2))^2+(2t\cos(t^2))^2}dt\\
= \int_0^\tau 2t \sqrt{\sin^2{(t^2)}+\cos^2{(t^2)}} dt\\
= \int_0^\tau 2t dt\\
= [t^2]_0^\tau\\
= \tau^2
$$

# Exercise 2

## Exercise 2-1

曲線 
$$
x = \cos{t}\\
y = \sin{t}\\
t \in [0,2\pi)
$$
を描き、その曲線の上に適当に多数の点を取り、その点での進行方向ベクトル$(\frac{dx}{dt},\frac{dy}{dt})$を描け。

Draw the curve and put its speed vectors $(\frac{dx}{dt},\frac{dy}{dt})$ on many points on the curve.

## Exercise 2-2

$$
x = \cos{t^2}\\
y = \sin{t^2}\\
t \in [0,\sqrt{2\pi})
$$
について同様のことをせよ。

Do the same for this parameterization.

## Exercise 2-3

3次元空間の曲線
$$
x = r\cos{t}\\
y = r\sin{t}\\
z = r^2\\
r = 0.1t+0.1\\
t \in [0,100]
$$
```{r,rgl=TRUE}
t <- seq(from=0,to=1,length=1000)*100
r <- 0.1*t+0.1
x <- r*cos(t)
y <- r*sin(t)
z <- r^2
plot3d(x,y,z,type="l")
```

$\frac{dx}{dt},\frac{dy}{dt},\frac{dz}{dt}$を計算し、曲線の長さとしてその積分をRのintegrate()関数を使って計算せよ。

This is a curve in 3D space. Calculate $\frac{dx}{dt},\frac{dy}{dt},\frac{dz}{dt}$ and calculate the length of curve with R's integrate() function.

また、曲線の長さを折れ線の長さの集まりとして
```{}
t <- seq(from=0,to=1,length=n)*100
```
$n=10,100,1000,10000,100000$について計算せよ。

For multiple n values, calculate the length of broken lines that approximate the curve.

# 単位円の面積 Area of a unit circle.

単位円の面積は、The following is the 1/4 of unit circle's area.

$$
y = \sqrt{(1-x^2)}\\
y = 0 \\
x \in [0,1]
$$
の4倍として求めることができる。

$$
4 \int_0^1 \sqrt{(1-x^2)}dx
$$

この積分を解くのは難しい。
It's no easy to solve this.

```{r}
integrate(function(x){4*sqrt(1-x^2)},0,1)
```

微小底辺$dl$を円周上に持ち、高さが半径1の三角形の和として計算することもできる。

中心からの半径方向を$\theta \in [0,2\pi)$とすれば

The area of circle is the sum of small triangles whose base lines are small arcs $dl$ of the circle and whose hight are the radius.
The directions of the triangles are parameterized as $\theta \in [0,2\pi)$.

$$
\int_0^{2\pi} \frac{1}{2}\times dl \\
dl = \sqrt{(\frac{dx}{d\theta})^2+(\frac{dy}{d\theta})^2}d\theta = d\theta 
$$
これは、円周の長さ$2\pi$を底辺の長さとする三角形の面積。

It corresponds to the area of triangle whose base is a whole circle.

別の方法は円を半径が0から1に変化する円のあつまりとみなすものである。

Another way is to consider the area of circle as the sum of cicles with radius ranging from 0 to 1.

半径$r$での微小面積は$2\pi r \times dr$.

The small area attached to the circle with radius r is $2\pi r \times dr$.

$$
\int_0^1 2\pi r dr = 2\pi \int_0^1 r dr = 2\pi [\frac{1}{2}r^2]^1_0 = \pi
$$
# 曲面 Curved surface

以下のような曲面の面積はどうなるか?

How do we approach to the area of the curved surface as below?

$$
z = \sin{x}+\cos{y}\\
x, y \in [-5,5]
$$
```{r,rgl=TRUE}
x <- y <- seq(from=-1,to=1,length=100) * 5
xy <- as.matrix(expand.grid(x,y))
z <- sin(xy[,1])+cos(xy[,2])
xlim <- max(xy[,1])-min(xy[,1])
ylim <- max(xy[,2])-min(xy[,2])
zlim <- max(z)-min(z)
plot3d(xy[,1],xy[,2],z,aspect=c(xlim,ylim,zlim))
```

(x,y),(x+dx,y),(x,y+dy),(x+dx,x+dy)が作る「斜面」の面積$dS$を足し合わせる。

Small areas that are slopes over (x,y),(x+dx,y),(x,y+dy),(x+dx,x+dy), $dS$, should be summed up.

これを斜めになった平行四辺形の面積とみなす。
They are considered small parallelograms.

平行四辺形の第1の辺は$(x,y,z)$と$(x+dx,y,z+\frac{\partial x}{\partial z}dx)$とが作るベクトル。
第2の辺は$(x,y,z)$と$(x,y+dy+z+\frac{\partial y}{\partial z}dy)$。

The parallelogram's two edges are; one connecting $(x,y,z)$と$(x+dx,y,z+\frac{\partial x}{\partial z}dx)$ and 
and the other connecting $(x,y,z)$ and $(x,y+dy+z+\frac{\partial y}{\partial z}dy)$.

2つのベクトルは They are:
$$
v_1 = (dx,0,\frac{\partial z}{\partial x}dx)\\
v_2 = (0,dy,\frac{\partial z}{\partial y}dy)
$$

2つのベクトル$v_1,v_2$が作る平行四辺形の面積は
The area of parallelogram is;

$$
dS = ||v_1||v_2|\sin{\theta}|
$$
$$
\cos{\theta} = \frac{(v_1,v_2)}{|v_1||v_2|}
$$
ただし、$(v_1,v_2)$は内積。
where $(v_1,v_2)$ is inner product.

したがって Therefore,

$$
dS = ||v_1||v_2|\sqrt{1-(\frac{(v_1,v_2)}{|v_1||v_2|})^2}|\\
= \sqrt{|v_1|^2|v_2|^2 - (v_1,v_2)^2}
$$

ここで Using 
$$
|v_1|^2 = (1+(\frac{\partial z}{\partial x})^2dx^2\\
|v_2|^2 = (1+(\frac{\partial z}{\partial y})^2dy^2\\
(v_1,v_2) =(\frac{\partial z}{\partial x}\frac{\partial z}{\partial y})^2dxdy
$$

以上より、Then,
$$
dS = \sqrt{1+(\frac{\partial z}{\partial x})^2+(\frac{\partial z}{\partial y})^2} dxdy
$$

## Exercise 3

### Exercise 3-1 

2つの式を比較せよ。

Compare two formulas.

$$
\int_{0}^3 \sqrt{dx^2+dy^2} = \int_0^3 \sqrt{(\frac{dx}{dx})^2+(\frac{dy}{dx})^2}dx = \sqrt{1+(\frac{dy}{dx})^2}dx \\
$$

$$
dS = \sqrt{1+(\frac{\partial z}{\partial x})^2+(\frac{\partial z}{\partial y})^2} dxdy
$$


## 曲面の前に、平面 Before using this to curved surface, apply it to the flat one.

$$
z = ax + by + c
$$

$$
\frac{\partial z}{\partial x}=a\\
\frac{\partial z}{\partial y}=b\\
dS = \sqrt{1+a^2+b^2}dxdy \\
\int_{x_0}^{x_1} \int_{y_0}^{y_1} dS dxdy = \sqrt{1+a^2+b^2} \times (x_1-x_0)(y_1-y_0)
$$
単なる、平行四辺形の面積。

Simple parallelogram's area.

## 単位球面積 Surface area of unit sphere.

$$
x^2+y^2+z^2=1\\
z = \pm \sqrt{1-(x^2+y^2)}
$$
$ z= \sqrt{1-(x^2+y^2)}$ , $0 \le x^2+y^2 \le 1;x,y \ge 0$を扱うことにする。

Calculate a part of the sphere of $ z= \sqrt{1-(x^2+y^2)}$ , $0 \le x^2+y^2 \le 1;x,y \ge 0$.

$$
dS = \sqrt{1+(\frac{\partial z}{\partial x})^2+(\frac{\partial z}{\partial y})^2} dxdy\\
\frac{\partial z}{\partial x}= \frac{1}{2}\frac{-2x}{\sqrt{1-(x^2+y^2)}}\\
\frac{\partial z}{\partial y} = \frac{1}{2}\frac{-2y}{\sqrt{1-(x^2+y^2)}}\\
$$

$$
dS = \sqrt{1+\frac{x^2}{1-(x^2+y^2)}+\frac{y^2}{1-(x^2+y^2)}}dxdy\\
= \sqrt{\frac{1}{1-(x^2+y^2)}}dxdy
$$

$$
S = \int_{x^2+y^2 \le 1} dS = \int_{x^2+y^2 \le 1} \sqrt{\frac{1}{1-(x^2+y^2)}}dxdy
$$

これは解けないが、うまくパラメタ表示を変えることで解けるようにできる。

This is not easy but converting to different parameterization, then it can  be.

$$
x = r \cos{\theta}\\
y = r \sin{\theta}\\
0 \le r \le 1\\
0 \le \theta \le 2\pi\\
$$

$(r\cos{\theta},r\sin{\theta}),((r+dr)\cos{\theta},(r+dr)\sin{\theta}),(r\cos{(\theta+d\theta)},r\sin{(\theta+d\theta)}),((r+dr)\cos{(\theta+d\theta)},(r+dr)\sin{(\theta+d\theta)})が作る微小面積は

Small area should be changed accordingly.

$$
dr \times (r \times d\theta)
$$
したがって$dxdy$に相当するのは$rdrd\theta$

Means $dxdy$ is converted to $rdrd\theta$.

ゆえに
$$
S = \int_{x^2+y^2 \le 1} \sqrt{\frac{1}{1-(x^2+y^2)}}dxdy\\
= \int_{0 \le r \le 1, 0\le \theta \le 2\pi} \sqrt{\frac{1}{1-r^2(\cos^2\theta + \sin^2{\theta})}}rdrd\theta\\
=\int_0^{2\pi} 1 d\theta \int_0^1 r\sqrt{\frac{1}{1-r^2}}dr
$$

$$
\frac{d}{dr} \sqrt{1-r^2} = \frac{1}{2}(-2r) \sqrt{\frac{1}{1-r^2}}\\
= -r\sqrt{\frac{1}{1-r^2}}
$$

したがって
$$
S = 2\pi
$$
これは球の上半分であるから、球面全体の表面積は、$4\pi$.

This is the upper half's area, therefore the are of unit sphere is $4\pi$.

## 曲面 Curverd surface


$$
z = \sin{x}+\cos{y}\\
x, y \in [-5,5]
$$
```{r,rgl=TRUE}
x <- y <- seq(from=-1,to=1,length=100) * 5
xy <- as.matrix(expand.grid(x,y))
z <- sin(xy[,1])+cos(xy[,2])
xlim <- max(xy[,1])-min(xy[,1])
ylim <- max(xy[,2])-min(xy[,2])
zlim <- max(z)-min(z)
plot3d(xy[,1],xy[,2],z,aspect=c(xlim,ylim,zlim))
```

$$
\frac{\partial z}{\partial x} = -\sin{x}\\
\frac{\partial z}{\partial y} = \cos{y}
$$


2つのベクトルは The vecors are;
$$
(dx,0,-\sin{x}dx)\\
(0,dy,\cos{y}dy)
$$


したがって Then

$$
dS = ||v_1||v_2|\sqrt{1-(\frac{(v_1,v_2)}{|v_1||v_2|})^2}|\\
= \sqrt{|v_1|^2|v_2|^2 - (v_1,v_2)^2}
$$

ここで Using
$$
|v_1|^2 = (1+\sin^2{x})dx^2\\
|v_2|^2 = (1+\cos^2{y})dy^2\\
(v_1,v_2) = -\sin{x}\cos{y}dxdy
$$

以上より、Therefere
$$
dS = \sqrt{(1+\sin^2{x})(1+\cos^2{y})-\sin^2{x}\cos^2{y}} dxdy
= \sqrt{1+(\sin^2{x}+\cos^2{y})}dxdy
$$

解けないけれど、計算機的には計算できる。
Computationally calculable.


## Exercise 4

## Exercise 4-1

Mathematicaのサイトを使ってこの面積を計算せよ。

Calculate the area using Mathematica's site.

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

[][][][]多変量ベイズ推定のための偏微分

---
title: "Calculus4 Partial differentiation"
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)
```

# 正規分布のパラメタ推定 Estimation of parameters of normal distribution

$k$個の実数 $X=(x_1,...,x_k)$が観察されたとする。

Assume k real values $X=(x_1,...,x_k)$, are observed.

これらが平均$m$、SD $s$の正規分布からの独立標本であるとすると、その尤度関数は

The likelihood function under the hypothesis where they are from normal dist with mean $m$ and SD $s$ is;

$$
L(m,s|x) = \prod_{i=1}^k (\frac{1}{\sqrt{2\pi s^2}}e^{-\frac{1}{2}\frac{(x_1-m)^2}{s^2}})
$$
である。

対数をとって、定数部分を省略すれば

Taking logaithm,

$$
\log{L} = -k\log{s} -\frac{1}{2}\frac{1}{s^2}\sum_{i=1}^k (x_i-m)^2 + C
$$

これを、$m$,$s$の関数らしく変形すると

Transform the function so that it appears a function of $m$ and $s$,

$$
\log{L} = -k\log{s}-\frac{1}{2}\frac{1}{s^2} (k m^2 - 2(\sum_{i=1}^k x_i)m + \sum_{i=1}^k x_i^2) + C
$$

この対数尤度に準じた関数を、横軸に$m$、縦軸に$s$をとって描いてみる。

Draw the function with $m$ horizontal and $s$ vertical axes;

```{r,rgl=TRUE}
k <- 10
m0 <- 2
s0 <- 1
x <- rnorm(k,m0,s0)

m <- seq(from=-2,to=4,by=0.05)
s <- seq(from=0.5,to=20,by=0.05)
ms <- as.matrix(expand.grid(m,s))
log.L <- -k*log(ms[,2])-1/2 * (1/ms[,2]^2) * (k*ms[,1]^2-2*sum(x)*ms[,1]+sum(x^2))
plot3d(ms[,1],ms[,2],log.L)
```

## m = 0

平均が0の正規分布からのサンプルであることがわかっているなら

Under the assumption where mean is 0,

```{r,rgl=TRUE}
id <- which(ms[,1]==0)
plot3d(ms[,1],ms[,2],log.L)
spheres3d(ms[id,1],ms[id,2],log.L[id],radius=3,color=2)
```

$m=0$の下での、$s$の最尤推定値は、対数尤度(に準じた)関数の$m$に0を代入した尤度関数を微分すればよい

The function with $m=0$ should be differentiate to tell the MLE of $s$. 

$$
\log{L} = -k\log{s}-\frac{1}{2}\frac{1}{s^2} (k m^2 - 2(\sum_{i=1}^k x_i)m + \sum_{i=1}^k x_i^2) + C\\
=-k\log{s}-\frac{1}{2}\frac{1}{s^2} \sum_{i=1}^k x_i^2 + C
$$
$$
\frac{d}{ds} (\log{L}) = -\frac{k}{s} + \frac{1}{s^3} \sum_{i=1}^k x_i^2\\
=-\frac{1}{s^3}(k s^2 - \sum_{i=1}^k x_i^2)
$$

結局 then,
$$
s = \sqrt{\frac{\sum_{i=1}^k x_i^2}{k}}
$$

```{r}
id <- which(ms[,1]==0)
plot3d(ms[,1],ms[,2],log.L)
spheres3d(ms[id,1],ms[id,2],log.L[id],radius=3,color=2)
s. <- sqrt(sum(x^2)/k)
log.L. <- -k*log(s.)-1/2 * (1/s.^2) * (sum(x^2))
spheres3d(0,s.,log.L.,radius=5,color=3)

```

## m=1

$m=1$で同じことをするなら

In the case of $m=1$,
$$
\log{L} = -k\log{s}-\frac{1}{2}\frac{1}{s^2} (k m^2 - 2(\sum_{i=1}^k x_i)m + \sum_{i=1}^k x_i^2) + C\\
=-k\log{s}-\frac{1}{2}\frac{1}{s^2} (k-2\sum_{i=1}^k x_i + \sum_{i=1}^k x_i^2) + C
$$

$$
\frac{d}{ds} (\log{L}) = -\frac{k}{s} + \frac{1}{s^3} (k-2\sum_{i=1}^k x_i + \sum_{i=1}^k x_i^2)\\
=-\frac{1}{s^3}(k s^2 - \sum_{i=1}^k (k-2\sum_{i=1}^k x_i + \sum_{i=1}^k x_i^2))
$$
$$
s = \sqrt{\frac{\sum_{i=1}^k (k-2\sum_{i=1}^k x_i + \sum_{i=1}^k x_i^2)}{k}}
$$
```{r}
id <- which(ms[,1]==0)
plot3d(ms[,1],ms[,2],log.L)
spheres3d(ms[id,1],ms[id,2],log.L[id],radius=3,color=2)
s. <- sqrt(sum(x^2)/k)
log.L. <- -k*log(s.)-1/2 * (1/s.^2) * (sum(x^2))
spheres3d(0,s.,log.L.,radius=5,color=3)

id <- which(ms[,1]==1)
spheres3d(ms[id,1],ms[id,2],log.L[id],radius=3,color=4)
s.. <- sqrt((k-2*sum(x)+sum(x^2))/k)
log.L.. <- -k*log(s..)-1/2 * (1/s..^2) * (k-2*sum(x)+sum(x^2))
spheres3d(1,s..,log.L..,radius=5,color=5)
```

## $m= ? $

$m=0$, $m=1$の場合では、$m$を固定して$\log{L}$を$s$で微分した。

For the cases of $m=0$ and $m=1$, $\log{L}$ was differentiated with $s$ with $m$ as constant.

このように変数を固定して残りの変数で微分するのが偏微分。
This is partial differentiation.

$$
\frac{\partial}{\partial s} \log{L} = -\frac{k}{s} + \frac{1}{s^3} (k m^2 - 2(\sum_{i=1}^k x_i)m + \sum_{i=1}^k x_i^2)
$$
$m$の値に応じて、$s$の最尤推定値は次のような値をとることがわかる

Partial differentiation tells MLE for arbitrary $m$.

$$
\sqrt{\frac{(k m^2 - 2(\sum_{i=1}^k x_i)m + \sum_{i=1}^k x_i^2)}{k}}
$$

```{r}
s.mle <- sqrt((k*m^2-2*sum(x)*m+sum(x^2))/k)

plot(m,s.mle,type="l")
```

## Exercise 1

### Exercise 1-1
$m$に関する偏微分をして、$s$とそれに対応する$m$の最尤推定値の関係をプロットせよ1

Partially differentiate $\log{L}$ with $m$ and draw the relation between $s$ and MLE of $m$.


### Exercise 1-3

常染色体上のSNPは3種類のディプロタイプを作る。
そのディプロタイプ頻度を$(p_{MM},p_{Mm},p_{mm})$とする。

SNP in autosomal chromosome makes 3 diplitypes, whose frequency is 
$(p_{MM},p_{Mm},p_{mm})$.


今、3ディプロタイプ人数が$(n_{MM},n_{Mm},n_{mm})$と観察されたとする。

Assume $(n_{MM},n_{Mm},n_{mm})$ are the observed number of individuals of three diplotypes.

この観察の下での、集団のディプロタイプ頻度
(p_{MM},p_{Mm},p_{mm})の尤度関数は

The likelihood fucntion of diplotype frequency for the observarion follows.

$$
L((p_{MM},p_{Mm},p_{mm})|(n_{MM},n_{Mm},n_{mm})) = \begin{pmatrix} n_{MM}+n_{Mm}+n_{mm}\\n_{MM},n_{Mm},n_{mm} \end{pmatrix} p_{MM}^{n_MM} p_{Mm}^{n_{Mm}} p_{mm}^{n_{mm}}
$$


尤度関数の対数を取り、定数部分を省略すると以下のようになる。

Taking logarithm,

$$
\log{L((p_{MM},p_{Mm},p_{mm})|(n_{MM},n_{Mm},n_{mm}))} = n_{MM} \log{p_{MM}} + n_{Mm}\log{p_{Mm}} + n_{mm}\log{p_{mm}} + C
$$

$p_{mm}=1 - p_{MM}-p_{Mm}$であることを利用して、$\log{L}$を$(p_{MM},p_{Mm})$の2変数関数であるとみなし、$p_{MM}$,$p_{Mm}$でそれぞれ偏微分せよ。

Using $p_{mm}=1 - p_{MM}-p_{Mm}$, handle the function as the funtion with two parameters $p_{MM}$,$p_{Mm}$ and partially differentialte the function with both.



### Exercise 1-4

Exercise 1-3 の偏微分を利用して、$(p_{MM},p_{Mm},p_{mm})$の最尤推定値を求めよ

Using the result of Exercise 1-4, answer MLE of $(p_{MM},p_{Mm},p_{mm})$.

### Exercise 1-5

アレル頻度を(p,1-p)とすると、Hardy-Weinberg 平衡の下でのディプロタイプ頻度は

Under the condition of Hardy-Weinberg equilibrium, three diplotype freq when allele frequencies are (p,1-p),

$$
(p_{MM},p_{Mm},p_{mm}) = (p^2, 2p(1-p),(1-p)^2)
$$
である。HWEを仮定すると、尤度関数は、1変数関数となる。

Under the assumption of HWE, the likelihood function is a funciton with one parameter.

HWE仮定の下での、尤度関数を示し、それを微分することでアレル頻度$p$の最尤推定値を求めよ。

Show the likelihood function and differentiate it and answer MLE.

### Exercise 1-6

Hardy-Weinberg不平衡を許せば

Under the assumption of HW-disequiriblium,


$$
(p_{MM} + \delta, p_{Mm} - 2\delta, p_{mm} + \delta)
$$
と表せる。

これにより、尤度関数は$p$,$\delta$の2変数関数となる

2変数で偏微分し、その結果を示せ

Now the likelihood funtion has two parameters.
Partially differentiate it and show the results.

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

[][][][][]確率密度分布、累積分布、生存分布、ハザード分布

---
title: "Calculus3 Density, Cumulative, Hazard"
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)
```


確率密度関数を積分すると累積密度関数。

Integrate probability density function, then it is cumulative density function.

累積密度関数を微分すると確率密度関数。
Differentiate cumulative density function, then it is density function.

確率分布の本体は同じで、その見方が異なっているだけ。
Both density and cumulative density functions represent the same function and they are two ways to show it.

いろいろな見方を扱い、それらの微積分関係を確認する。

Some other ways to represent probability function follow.

# 確率密度関数 Probability density function $f(x)$

確率変数$X$が値$x$をとる確率密度を$P(X=x)$と書く。

Let $P(X=x)$ denote the probability density of X being x 

確率密度関数 Probability density function.
$$
f(x) = P(X=x)\\
f(x) \ge 0\\
\int_{x\in D} f(x) = 1
$$

指数分布を例に取ると In case of exponential distribution,

$$
f(x) = \lambda e^{-\lambda x}
D = [0;\infty)
$$

```{r}
lambda <- 3
x <- seq(from=0,to=5,length=100)
fx <- lambda * exp(-lambda *x)
plot(x,fx,type="l")
fx. <- dexp(x,lambda)
plot(fx,fx.)
```

# 累積分布関数 Cumulative distribution function $F(x)$

$$
F(x) = P(X \le x) = \int_{}^x f(x) dx\\
\frac{d}{dx} F(x) = f(x)
$$

指数分布の場合 Exponential distribution:

$$
F(x) = \int_{0}^x f(t) dt\\
= \int_0^x \lambda e^{-\lambda t}dt\\
= [\lambda \frac{1}{-\lambda} e^{-\lambda t}]^x_0 \\
= - e^{-\lambda x} + e^{0}\\
= 1-e^{-\lambda x} 
$$
```{r}
Fx <- 1-exp(-lambda*x)
plot(x,Fx,type="l")
abline(h=1,col=2)
Fx. <- pexp(x,lambda,lower.tail=TRUE)
plot(Fx,Fx.)
```

# 生存関数 Survival function $S(x)$

The order of x in $F(x)$ is inversed.

累積分布関数$F(x)$の昇降順を逆にした関数。

$P(X=x)$を時刻$X=x$での瞬間死亡率とすると、$S(x)$は時刻$X=x$での生存者の割合

When the $P(X=x)$ is probability of death at time $X=x$, then $S(x)$ stands for the fraction of survivers at time $X=x$.

$$
S(x) = P(X>x)\\
S(x) = 1- F(x) = \int_x^{} f(x) dx\\
=1- (1-e^{-\lambda x}) \\
= e^{-\lambda x}
$$
```{r}
Sx <- exp(-lambda * x)
plot(x,Sx,type="l")
Sx. <- pexp(x,lambda,lower.tail=FALSE)
plot(Sx,Sx.)
```
# 逆累積分布関数 Inverse of cumulative density function $G(t)$



$$
P(X \le G(t)) = t\\
\int_{}^{G(t)} f(x) dx = t\\
F(G(t)) = t\\
G(t) = F^{-1}(t)
$$

$$
F(G(t)) = t\\
F(G(t)) = t = 1-e^{-\lambda G(t)}\\
e^{-\lambda G(t)} = 1-t\\
-\lambda G(t) = \log{1-t}\\
G(t) = -\frac{\log{(1-t)}}{\lambda}
$$


```{r}
t <- seq(from=0,to=1,length=100)
Gt <- -(log(1-t))/lambda
plot(t,Gt,type="l")
Gt. <- qexp(t,lambda)
plot(Gt,Gt.)
```

# 逆生存関数 Inverse of survival function $Z(x)$

$$
P(X>Z(t)) = t\\
\int_{Z(t)}^{} f(x) dx = t\\
S(Z(t)) = t\\
Z(t) = S^{-1}(t)
$$

$$
S(Z(t)) = t\\
S(Z(t)) = t = e^{-\lambda Z(t)}\\
\lambda Z(t) = - \log{t}\\
Z(t) = -\frac{\log{t}}{\lambda}
$$

```{r}
Zt <- -log(t)/lambda
plot(t,Zt,type="l")
Zt. <- qexp(t,lambda,lower.tail=FALSE)
plot(Zt,Zt.)
```

# ハザード関数 Hazard function $h(x)$

$$
h(x) = \frac{f(x)}{S(x)}\\
h(x) = \frac{f(x)}{\int_x^{} f(x) dx}
$$
$$
\frac{d}{dx} F(x) = f(x)\\
S(x) = 1- F(x)\\
\frac{d}{dx} S(x) = -\frac{d}{dx}F(x) = -f(x)
$$

$$
h(x) = -\frac{\frac{d}{dx}S(x)}{S(x)}\\
= - \frac{d}{dx} (\log{S(x)})
$$

$$
f(x) = \lambda e^{-\lambda x}\\
S(x) = e^{-\lambda  x}\\
h(x) = \frac{f(x)}{S(x)} = \lambda\\
h(x) = -\frac{d}{dx} \log{S(x)}\\
= -\frac{d}{dx} \log{e^{-\lambda x}}\\
= -\frac{d}{dx} (-\lambda x)\\
= \lambda
$$
生き残っている部分$S(x)$の死亡率(瞬間死亡率)

Probability to die for the survivers.

```{r}
hx <- fx/Sx
plot(x,hx,type="l")
```



指数分布とは、「どの時点においても、生き残っている部分での瞬間死亡率が一定」である分布

Exponential distribution stands for the stochastic process where the probability to die among the people alive is constant regardless of the size of population.

# 累積ハザード関数 Cumulative hazard function $H(x)$

$$
\frac{d}{dx}H(x) = h(x) = -\frac{d}{dx} \log{S(x)}\\
H(x) = -\log{S(x)}
$$

$$
H(x) = -\log{S(x)} = -\log{e^{-\lambda x}}= \lambda x
$$

# Exercise 1

## Exercise 1-1

生存時間のモデル、機械が故障するまでの時間のモデルのひとつが指数分布である。

Exponential distribution is one of stochastic models of survival time (time to death) or time to failure.

別のモデルとしてワイブル分布がある。

Weibull distribution is one of them.

ワイブル分布の確率密度分布は The probability density distribution of Weibull distribution is;

$$
f(x) = k\lambda (\lambda x)^{k-1} e^{-(\lambda x)^k}
$$
で与えられる。

このワイブル分布の生存関数、ハザード関数を求め、指数分布のそれと比較せよ

Calculate survival function and hazard function for Weibull distribution, then compare them with ones of exponential distribution.