ニュートン法をPythonで実装【統計学（最小二乗法）への応用】

## ニュートン法で学ぶ！方程式の解き方と統計学への応用
この記事では、ニュートン法の基礎を紹介します。
統計学への応用を視野に、統計学での例(最小二乗法)を紹介します。
しくみを実感できるように、Pythonのコードも紹介します。


## ニュートン法とは？
**ニュートン法**は、方程式 $ f(x) = 0 $ の解を効率的に見つける「数値計算」の方法です。

例：$ x^2 = 2 $ の解$ \sqrt{2} $を数値計算で近似するときに使えます。


---

## ニュートン法の納得できるしくみ
### 3ステップで理解！
1. **初期値設定**：適当な値 $ x_0 $ を選ぶ（例：$ x_0 = 1 $）
2. **接線を引く**：$ f(x) = x^2 - 2 $ の $ x_0 $ での接線を描く
3. **交点を更新**：接線とx軸の交点を次の近似値 $ x_1 $ とする

![](https://upload.wikimedia.org/wikipedia/commons/e/e0/NewtonIteration_Ani.gif)
（※図解：接線を繰り返し引いて解に近づく様子）
Wikiからの[引用](https://upload.wikimedia.org/wikipedia/commons/e/e0/NewtonIteration_Ani.gif)です。

### 計算式
$$ 
x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}
$$ 
- $ f'(x) $ は「関数の傾き」（例：$ f(x)=x^2-2 $ なら $ f'(x)=2x $）

### 具体例（√2の計算）
- $ x_0 = 1 $ → $ x_1 = 1.5 $ → $ x_2 ≈ 1.4167 $ → 3回で約1.4142（$ \sqrt{2} $）に！

---

## なぜ理論計算ではなくニュートン法を使うのか？
### 主な理由
1. **複雑な方程式**：$ x^5 + 3x - 2 = 0 $ など理論解が存在しない場合
2. **高速計算**：10回程度の計算で高精度（例：√2を小数第10位まで計算可能）


### 理論計算 vs ニュートン法
| 方法 | メリット | デメリット |
|------|----------|------------|
| 理論解（例：二次方程式の解の公式） | 正確な解が得られる | 複雑な方程式では不可能 |
| ニュートン法 | ほとんどの方程式に適用可能 | 近似値しか得られない |

---

## 最小二乗法とは？
最小二乗法は、データに最もよく当てはまる直線や曲線を見つける方法です。  
例えば「勉強時間とテストの点数の関係」を表す直線 $y = ax + b$ の傾き $a$ と切片 $b$ を求めるために使われます。  
**目標**：データとの誤差（予測値と実際の値の差）を最小化するパラメータを見つけることです。

---

## ニュートン法で最小二乗法を解く理由
最小二乗法は通常、数学的に解くことができますが、ニュートン法を使うと以下のメリットがあります：
1. **複雑な問題に対応可能**：非線形モデルや多変数モデルでも適用できる。
2. **高速収束**：少ない計算回数で高精度な解が得られる。

---

## 例題
### 問題
以下のデータに最もよく当てはまる直線 $y = ax + b$ の傾き $a$ と切片 $b$ を求めます。

| 勉強時間（時間） | テストの点数 |
|----------------|------------|
| 1              | 50         |
| 2              | 60         |
| 3              | 70         |

### 目標
誤差の二乗和（$E = \sum_{i=1}^n (y_i - (a x_i + b))^2$）を最小化する $a$ と $b$ を見つける。

---

## ニュートン法の適用
**誤差関数**：$E(a, b) = \sum_{i=1}^n (y_i - (a x_i + b))^2$

**勾配（一次導関数）**：
   - $\frac{\partial E}{\partial a} = -2 \sum_{i=1}^n x_i (y_i - (a x_i + b))$
   - $\frac{\partial E}{\partial b} = -2 \sum_{i=1}^n (y_i - (a x_i + b))$

**ヘッセ行列（二次導関数）**：
   - $\frac{\partial^2 E}{\partial a^2} = 2 \sum_{i=1}^n x_i^2$
   - $\frac{\partial^2 E}{\partial b^2} = 2n$
   - $\frac{\partial^2 E}{\partial a \partial b} = 2 \sum_{i=1}^n x_i$

ニュートン法では、これらの勾配とヘッセ行列を使ってパラメータ $a$ と $b$ を更新します。

更新式
```math
\begin{bmatrix}
a_{\text{new}} \\
b_{\text{new}}
\end{bmatrix}
= \begin{bmatrix}
a \\
b
\end{bmatrix}
- H^{-1} \nabla E
```

---

## Pythonで解く
以下は、ニュートン法を使って最小二乗法を解くPythonコードです。

```python
import numpy as np

# データ
X = np.array([1, 2, 3])  # 勉強時間
y = np.array([50, 60, 70])  # テストの点数
n = len(X)

# 初期値
a, b = 0.0, 0.0

# ニュートン法の反復
for _ in range(100):
    # 予測値
    y_pred = a * X + b
    
    # 勾配（一次導関数）
    grad_a = -2 * np.sum(X * (y - y_pred))
    grad_b = -2 * np.sum(y - y_pred)
    
    # ヘッセ行列（二次導関数）
    hess_aa = 2 * np.sum(X**2)
    hess_bb = 2 * n
    hess_ab = 2 * np.sum(X)
    
    # ヘッセ行列の逆行列
    det = hess_aa * hess_bb - hess_ab**2
    inv_hess_aa = hess_bb / det
    inv_hess_bb = hess_aa / det
    inv_hess_ab = -hess_ab / det
    
    # パラメータ更新
    delta_a = -(inv_hess_aa * grad_a + inv_hess_ab * grad_b)
    delta_b = -(inv_hess_ab * grad_a + inv_hess_bb * grad_b)
    a += delta_a
    b += delta_b

print(f"最適な直線: y = {a:.1f}x + {b:.1f}")
# 出力: 最適な直線: y = 10.0x + 40.0
```

## 結果の解釈
傾き $ a=10.0 $
勉強時間が1時間増えると、テストの点数が10点上がる。

切片 $ b=40.0 $
勉強時間が0時間の場合、テストの点数は40点と予測される。

この直線$ y=10x+40 $は、データに最もよく当てはまる直線であると言えます。