ぱたへね

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

F#でホデリングT^2法

入門機械学習による異常検知の最初の例題をF#で書いてみました。

www.amazon.co.jp

内容は入力されたデータの分布推定を行い、異常度の計算を行います。詳細は本を買ってください。

異常度の式 a(x') = ( \frac {x' - \hat \mu }{ \hat \sigma} ) ^ 2 を使ってデータから異常検知を行います。

F#での実行結果。Rと同じ結果になりました。

f:id:natsutan:20200807122957p:plain
異常検知

詰まりそうな所は2箇所です。

Rのライブラリをcsv化

Rを立ち上げて、write.csv()を使えばOK

> library(car)
> data(Davis)
> write.csv(Davis, "C:/home/tmp/car.csv")

カイの二乗分布の計算

次にRのqchisqに相当する関数をどうやって用意するかですが、MathnetのChiSquaredを使いました。

numerics.mathdotnet.com

qchisqに相当する関数はInvCDFなので引数に注意して呼び出せばOKです。

ソースコード

ソースはこんな感じで。

// Learn more about F# at https://fsharp.org
// See the 'F# Tutorial' project for more help.
open FSharp.Data
open FSharp.Charting
open MathNet.Numerics.Distributions

let csv_path = @"C:\home\myproj\study_ml\abnormaly_detection\data\car.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/car.csv">

let getWeights = 
    let car = CarDB.Load(csv_path)
    
    [ for row in car.Rows -> (float row.Weight)]
    
let plotWeights (xs :List<float>)  = 
    let n = Seq.length xs
    let chart =  Chart.Point [for i in 0 .. n - 1 -> (i + 1,  xs.[i]) ]
    let output_file = output_path + @"weights.png" 
    Chart.Save output_file chart
    
    printf "%s\n"  output_file
    None

let statWeight (xs:List<float>) =
    let mu = List.average xs
    let sigma2 = List.map (fun x -> (x - mu) ** 2.0) xs  |> List.average

    (mu, sigma2)

let calcAbnormality (xs:List<float>) =
    let mu, sigma2 = statWeight xs
    List.map (fun x -> (x - mu) ** 2.0 / sigma2) xs
   
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


[<EntryPoint>]
let main argv =
    let weights = getWeights
    plotWeights weights |> ignore

    let abnormality = calcAbnormality weights
    let threshhold =  ChiSquared.InvCDF(1.0, 0.99)
    printfn "閾値 = %A" threshhold

    plotAbnormality (abnormality, threshhold) |> ignore


    0 // return an integer exit code