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

---
title: "Calculus6 Taylor expansion テイラー展開"
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-値を得るとする。

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

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 を読み、積率母関数とテイラー展開との関係を説明せよ。

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

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)


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)



}
`