データサイエンスと生命科学

ニュートン法とは、実数値関数の根を数値的に計算する代表的なアルゴリズムです。

実数値関数の根を求めるとは、例えば以下のような方程式の解を求めることです。

$$ x^{3} + x - 8 = 0 $$

関数の根とは、f(x) = 0になるxの値のことです。たとえば、中学校の数学の方程式の授業では、x² - 4 = 0を解け、というような問題が出ます。答えは、x = 2またはx = -2のときです。これらの点が、関数f(x) = x² - 4の根になります。

方程式の解き方には、解析的に解くやり方と数値的に解くやり方があります。解析的に解く方法は、数式をきちんと立て、その数式を変形して正確な解を示すことです。たとえば、関数f(x) = x² - 4の場合、式変形によって正確な解を求めることができるので、解析的な計算ができます。中学校や高校で習う方程式の計算は、基本的には全てこの解析的な解を求めています。でも、方程式が複雑な場合は、解析的に解くことが難しい場合もあります。

数値的に解くとは、方程式の解を具体的な数字として計算機を使って力づくで答えを導き出す方法です。例えば、いろいろな値をxに代入してみて、しらみつぶしに探す、というのは最もシンプルな数値的な解き方ですが、効率が悪いです。そのため、情報科学では、効率的に計算するためのアルゴリズムが研究されており、ニュートン法はその１つです。数値的に解くときには、演算回数が少なく精度がよいことなどが望ましいです。

ニュートン法では、初めに適当な推測値からスタートして、より正確な解の近似値を探索します。反復を繰り返すことにより、正解に非常に近い値を得ることができます。

例として、以下の３次方程式の解をニュートン法を使ってどのように解くか考えてみましょう。

$$ x^{3} + x - 8 = 0 $$

この関数 (f(x) = x³ + x - 8)をRで視覚化してみると、方程式の解であるf(x) =0となるxは2の周辺にあるように見えます。

#Rでf(x)を描画してみる
x <- seq(-3, 3, length=100)
fx <- x^3 + x - 8
plot(x, fx, type="l")
abline(h=0,col="red")

f(x) = x³ + x - 8を描画した図。赤線はf(x) = 0を表す。交点のx座標の値から、f(x)=0を満たすxの値は2の近くであることがわかる。

ニュートン法は必ず解が求める保証はありませんが（その条件は複雑）、真の解に近い初期値を使えば多くの場合収束します。上の描画より、初期値x_0を2としてニュートン法を使って解を探します。

初期値が決まったら、反復プロセスを使用してより良い正確な値を計算していきます。ニュートン法では、$x_n$が現在の暫定解である場合、次の暫定解$x_{n+1}$は以下の式によって与えられます。

$$ x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} $$

f'(x)は、対象となる関数を微分して得られた導関数で、以下のように求められます。

$$ f(x) = x^{3} + x - 8 \\ f'(x) = 3x^{2} + 1 \\ $$

現在、初期値x_0=2が与えれれているので、上の式を使ってx_1が計算できます。

$$ x_{1} = x_0 - \frac{f(x_0)}{f'(x_0)} $$

同様に、x_2は

$$ x_{2} = x_1- \frac{f(x_1)}{f'(x_1)} $$

です。

これを繰り返していくと、どんどん真の解に近づいていく、ということです。そのため、初期値x_0, f(x) , f'(x)がきまれば、コンピューターで計算を繰り返していけば解を探索できます。本当にできるか、実際にRでやってみます。

f <- function(x){x^3 + x - 8}
fd <- function(x){3 * x^2 + 1} #f'(x)
n_iter <- 10
xn <- 2 #初期値x_0を設定する
xn_hist <- NULL
for(i in 1:n_iter){
    xn_hist <- c(xn_hist, xn)
    xn <- xn - (f(xn)/fd(xn))
}
plot(xn_hist, xlab="n", ylab="x_n", cex.lab=2)

x_nの値の変化。x_0は初期値の2であるが、nが増えていくごとに一定の値に収束していく。

x_nのnの値が大きくなるごとに、一定の値に収束していきます。 10回反復すると、x_10の値は1.833751になっていました。この時のf(x)はゼロになっていることも以下のようなコマンドで確かめられます。

#Solution
cat("x_n: ", xn, "\n")
cat("f(x): ", f(xn))

以上より、方程式x³ + x - 8 = 0の解は、x = 1.833751と計算できました。

RでBioconductorのmethylclockを使いたい。

www.bioconductor.org

まず、Docker上で最新版のR環境を構築。

hub.docker.com

docker run -it --name "methylclock" r-base /bin/bash

でRを起動する。

しかし、ここでいきなりインストールしようとしても、足りないパッケージが多くてうまくいかない。そこでまず、以下をインストールする。

apt install libxml2-dev

apt install libcurl4-openssl-dev

apt install cmake

apt install libssl-dev

apt install libfontconfig1-dev

apt install libharfbuzz-dev libfribidi-dev

apt install libfreetype6-dev libpng-dev libtiff5-dev libjpeg-dev

そのあと、R上で以下コマンドでインストールできる。

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("methylclock")

データサイエンスと生命科学

ニュートン法を使って方程式を解く[R][データサイエンス]

DockerとRでmethylclock packageをインストール

R3.6にscDDパッケージをインストール

LinuxでのTeamviewerの自動起動

debファイルからのインストールはgdebiを使う。

arXiv備忘録

任意のサイズで画像切り出し