F#で多変量正規分布に基づく異常検知

natsutan.hatenablog.com

の続き。入門機械学習による異常検知の多変量正規分布に基づく異常検知をF#で書いてみました。

入門機械学習による異常検知―Rによる実践ガイド | 剛, 井手 |本 | 通販 | Amazon

出てくる数式とキーワードはこんな感じで、詳細は本を買って読んでみてください。

平均値

$\hat{{\bf { \mu }}} = \frac 1 N \sum_{n=1}^N {\bf { x }}^{(n)}$

共分散行列

$\Sigma = \hat{\Sigma} = \frac 1 N \sum_{n=1}^N ({\bf { x }}^{(n)} - {\bf {\mu}}) ({\bf { x }}^{(n)} -\hat{\bf { \mu }}) ^T$

異常度の定義、マハラノビス距離

$a({\bf { x }}') = ({\bf { x }}' - \hat { \mu } ) ^ T \hat { \Sigma } ^ {-1}({\bf { x }} ' - \hat{{\bf { \mu }}})$

実行結果

教科書と同じ結果になりました。

f:id:natsutan:20200809164626p:plain

F#実装

MathNet周りでちょっと詰まったところを整理しました。

List、配列からMatrix型を作る

一度array2Dを作って、Matrix.Build.DenseOfArrayを呼び出せばOK。list0とlist1に同じ数のデータが入っているとして、このように書けばMatrixが作れます。

let arr = array2D [ list0; list1 ]
let data = Matrix.Build.DenseOfArray  arr

Matrix型からRowやColumnを取り出す。

Column(i), Row(i)を使う。2次元のMatrix xから0番目のコラムと1番目のコラムを取り出すにはこのようにします。

let c0 = x.Column(0)
let c1 = x.Column(1)

Matrixの平均値を求める。

MathNet.Numerics.Statistics.meanを使う。体重の値がはいったMatrixから平均値を求めるにはこうします。

let weight_mean = Statistics.Mean weights

RowやColumnでイテレートする

2次元MatrixからRow方向にデータを取り出して処理するには、EnumerateRowsを使います。

for row_data in  x.EnumerateRows() do
  ...

分からないこと

2次元配列のRow方向に一定の値を引く方法が良く分からないです。数式だと $x - \mu$ 　のように、体重や身長からそれぞれの平均値を引く操作です。

Pythonだと何も考えずにできます。

>>> import numpy as np
>>> a = np.array([[1,2,3], [4,5,6]])
>>> b = [[1],[2]]
>>> a - b
array([[0, 1, 2],
       [2, 3, 4]])
>>>

row方向にmapができれば良さそうだが、結局やり方分からず。

F#ソース

以下ソースです。

open FSharp.Data
open FSharp.Charting
open MathNet.Numerics.LinearAlgebra
open MathNet.Numerics.Statistics
open MathNet.Numerics.Distributions

let csv_path = @"C:\home\myproj\study_ml\abnormaly_detection\data\car2.csv"
let output_path = @"C:\home\myproj\study_ml\abnormaly_detection\fsharp\Hotelling\output\"
type CarDB = CsvProvider<"C:/home/myproj/study_ml/abnormaly_detection/data/car2.csv">

let plotAbnormality (xs :List<float>, th:float)  = 
    let mutable chart_list = []
    let n = Seq.length xs

    let chart1 =  Chart.Point [for i in 0 .. n - 1 -> (i + 1,  xs.[i]) ]
    let output_file = output_path + @"abnomality.png" 

    let chart2 = Chart.Line [ for x in 0.0 .. 200.0 -> (x, th) ]

    let chart_comb = Chart.Combine [chart1; chart2]

    Chart.Save output_file chart_comb
    
    printf "%s\n"  output_file
    None

let getData = 
    let car = CarDB.Load(csv_path)
        
    let weight = [ for row in car.Rows -> (float row.Weight)]
    let height = [ for row in car.Rows -> (float row.Height)]
    //weight, height
    let n = List.length weight

    //let data = DenseMatrix.init 2 n (fun i j -> if i == 0 then weight.Item j else height.Item j)
    let arr = array2D [ weight; height ]
    let data = Matrix.Build.DenseOfArray  arr

    data.Transpose()

let calcMahalanobis(x:Vector<float>, mu_hat:Vector<float>, sig_hat_inv:Matrix<float>) =
    let x_minus_mu_hat = x - mu_hat
    let a = x_minus_mu_hat * sig_hat_inv * x_minus_mu_hat

    a

let calcAbnormality (x:Matrix<float>) =
    let n = float x.RowCount
    let weights = x.Column(0)
    let heights = x.Column(1)

    // \hat μ
    let weight_mean = Statistics.Mean weights
    let height_mean = Statistics.Mean heights
    let mu_arr = [|weight_mean; height_mean|]
    let mu_hat = Vector.Build.DenseOfArray mu_arr

    // \hat Σ
    let sig_arr = array2D [ weights - weight_mean; heights - height_mean ]
    let x_minus_mu = Matrix.Build.DenseOfArray  sig_arr
    let sig_hat = (x_minus_mu * x_minus_mu.Transpose()) / n
    let sig_hat_inv = sig_hat.Inverse()

    // 異常度の計算
    let mutable abnormality = []
    for d in  x.EnumerateRows() do
        let a = calcMahalanobis(d, mu_hat, sig_hat_inv)
        abnormality <- a :: abnormality

    List.rev abnormality


[<EntryPoint>]
let main argv =    
    let data = hdtwo.getData
    let abnormality = hdtwo.calcAbnormality data
    
    let threshhold = ChiSquared.InvCDF(1.0, 0.99)
    hdtwo.plotAbnormality (abnormality, threshhold) |> ignore

    0 // return an integer exit code