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

---
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)
```

```
}