ぱたへね

はてなダイアリーはrustの色分けができないのでこっちに来た

F#正規方程式をで解く（行列計算）

F# 統計

統計的機械学習の数理100問 with Pythonに載っていた正規方程式をF#で解いてみました。F#で行列計算をするサンプルとして。

www.kyoritsu-pub.co.jp

正規方程式はこういう式です。

\hat{\beta} = (X^ TX)^ {-1}X^ Ty

Python 実装

教科書にのっているPython実装です。

import numpy as np

n = 100
p = 2
beta = np.array([1, 2, 3])
x = np.random.randn(n, 2)

y = beta[0] + beta[1] * x[:, 0] + beta[2] * x[:, 1] + np.random.randn(n)

X = np.insert(x, 0, 1, axis=1) # 左にすべて１の列を置く
beta_hat = np.linalg.inv(X.T @ X) @ X.T @ y

print(beta_hat)

実行するとこのように表示され、β=(1,2,3)が求められているようです。詳細知りたい人は本買ってください。

[0.91691516 1.91218893 2.91793258]

F#実装

open MathNet.Numerics.Distributions
open MathNet.Numerics.LinearAlgebra

let randn():double = 
    Normal(mean=0.0, stddev=1.0).Sample()
 
let makedata (N:int) = 
     // 一番左だけ1.0、残りはランダム
     let X = DenseMatrix.init 3 N (fun i j -> if i = 0 then 1.0 else randn())
     let beta = vector  [1.0; 2.0; 3.0]
     let noise = Vector<double>.Build.Random(N)
 
     let y = X.Transpose() * beta + noise
 
     X.Transpose(), y
 
let lsm (X:Matrix<double>, y:Vector<double>) =
     //正規方程式
     let beta = (X.Transpose() * X).Inverse() * (X.Transpose() * y)
 
     beta

実行するとだいたい合ってる感じの数値が出てきます。

beta = seq [1.002078054; 2.040203715; 3.098223575]

行列の操作、Row、Column単位のイテレート、map等、出来そうなんだけどやり方が分からないことがいっぱいなので、ぼちぼち勉強していきたい。