0

Baum-Welch アルゴリズムを実装し、既知の分布で生成されたおもちゃのデータで遊んでいます。データは正規分布しており、隠れた状態に応じて平均と標準偏差が異なります。2つの状態があります。アルゴリズムは、ランダム データに応じて常に (0; 1) または (1; 0) のいずれかに収束する隠れ状態の初期分布を除いて、ほとんどのパラメーターに対して収束するようです。

このアルゴリズムでは正常ですか？もしそうなら、バグを見つける方法のヒントではないにしても、いくつかの参考文献をいただければ幸いです。

コード (F#)。最初のヘルパー モジュール:

module MyMath

let sqr (x:float) = x*x

let inline (./) (array:float[]) (d:float) =
  Array.map (fun x -> x/d) array

let inline (.*) (array:float[]) (d:float) =
  Array.map (fun x -> x*d) array

let map f s =
  s |> Seq.map f |> Seq.toArray

let normalize v = 
  let sum = Seq.sum v
  map (fun x -> x/sum) v

let row i array = seq { for j in 0 .. (Array2D.length2 array)-1 do yield array.[i,j]}
let column j array = seq { for i in 0 .. (Array2D.length1 array)-1 do yield array.[i,j]}

let sum (v:float[]) = v |> Array.sum
let sumTo N (f:int->float) = Seq.init N f |> Seq.sum
let sum_column j (array:float[,]) = column j array |> Seq.sum
let sum_row i (array:float[,]) = row i array |> Seq.sum

let mean data = (sum data)/(float (Array.length data))
let var data = 
    let m=mean data
    let N=Array.length data
    let sum=Seq.sumBy (fun x -> sqr(x)) data
    sum/(float N)

let induction start T nextRow =
  let result =  Array.zeroCreate T
  result.[0] <- start
  for t=1 to T-1 do
    result.[t] <- nextRow t result.[t-1]
  result

let backInduction last T previousRow =
  let result =  Array.zeroCreate T
  result.[T-1] <- last
  for t=T-2 downto 0 do
    result.[t] <- previousRow t result.[t+1]
  result

let inductionNormalized start T nextRow =
  let result =  Array.zeroCreate T
  let norm = Array.zeroCreate T
  norm.[0] <- sum start
  result.[0] <- start./norm.[0]
  for t=1 to T-1 do
    result.[t] <- nextRow t result.[t-1]
    norm.[t] <- sum result.[t]
    result.[t] <- result.[t]./norm.[t]
  (result, norm)

メインモジュール:

module BaumWelch

open System
open MyMath

let mu (theta : float[,]) q = theta.[q,0]
let sigma (theta : float[,]) q = theta.[q,1]

let likelihood getDrift getVol dt parameters state observation  =
    let mu = getDrift parameters state
    let sigma =  Math.Abs (getVol parameters state:float)
    let sqrt_dt = Math.Sqrt dt
    let residueSquared = 
        let r = Likelihood.normalizedResidue mu sigma dt sqrt_dt observation in r*r
    let result = (Math.Exp (-0.5*residueSquared))/(sigma * (Math.Sqrt (2.0*Math.PI*dt)))
    if result<0.0 then failwith "Negative density, it certainly shouldn't have happened"
    else result

let alphaBeta b (initialPi:float[]) initialA observations= //notation in comments from the Erratum for Rabiner
    let T = Array.length observations
    let N = Array2D.length1 initialA
    let alphaStart = Array.init N (fun i -> initialPi.[i] * (b i observations.[0]))  //this contains \bar{\alpha}
    let alpha_j_t (previousRow:float[]) t j  = (sumTo N (fun i -> previousRow.[i]*initialA.[i, j]))* (b j observations.[t]) //this contains \bar{\alpha}
    let alphaInductionStep t previousRow = Array.init N (alpha_j_t previousRow t)
    let (alpha, norm) = inductionNormalized alphaStart T alphaInductionStep
    let betaStart = Array.init N (fun i -> 1.0/norm.[T-1])
    let beta_j_t (nextRow:float[]) t j = (sumTo N (fun i -> initialA.[j, i]*nextRow.[i]*(b i observations.[t+1])))/norm.[t]
    let betaInductionStep t nextRow = Array.init N (beta_j_t nextRow t)
    let beta = backInduction betaStart T betaInductionStep
    (alpha, beta, norm) //c_t = 1/norm_t


let log_P_O norm = 
  let result = norm |> Seq.sumBy (fun norm_t -> Math.Log norm_t)//c_t = 1/norm_t
  if Double.IsNaN result then failwith "log likelihood is NaN"
  else result

let gamma (alpha:float[][], beta:float[][], norm:float[])  i t = 
    alpha.[t].[i]*beta.[t].[i]*norm.[t]

let xi b (initialA:float[,]) (alpha:float[][]) (beta:float[][]) (observations:float[]) i j t = 
    alpha.[t].[i]*initialA.[i,j]*(b j observations.[t+1])*beta.[t+1].[j]


let oneStep llFunction dt (initialPi, initialA, initialTheta) observations =
  let T = Array.length observations
  let N = Array2D.length1 initialA
  let b = llFunction dt initialTheta
  let (alpha, beta, norm) = alphaBeta b initialPi initialA observations
  let gamma = gamma (alpha, beta, norm)
  let xi = xi b initialA alpha beta observations
  let pi = Array.init N (fun i -> gamma i 0) //Rabiner (40a)
  let A =  //Rabiner (40b)
    let A_func i j = (sumTo (T-1) (xi i j))/(sumTo (T-1) (gamma i))
    Array2D.init N N A_func
  let mean i = (sumTo T (fun t -> (gamma i t) * observations.[t]))/(sumTo T (gamma i))//Rabiner (53)
  let var i = 
    let numerator = sumTo T (fun t -> (gamma i t) * (sqr (observations.[t]-(mean i))))
    let denumerator = sumTo T (gamma i)
    numerator/denumerator
  let mu i = ((mean i) + 0.5*(var i))/dt
  let sigma i = Math.Sqrt ((var i)/dt)
  let theta = Array2D.init N 2 (fun i k -> if k=0 then mu i else sigma i) 
  let logLikelihood = log_P_O norm //Rabiner (103)
  (logLikelihood, (pi, A, theta))

let print (ll, (pi, A, theta)) = 
  printfn "pi = %A" pi
  printfn "A = %A" A
  printfn "theta = %A" theta
  printfn "logLikelihood = %f" ll

let baumWelch likelihood dt initialParams observations =
  let tolerance = 10e-5
  let rec doStep parameters previousLL =
    //print (previousLL, parameters)
    let (logLikelihood, parameters) = oneStep likelihood dt parameters observations
    if Math.Abs(previousLL - logLikelihood) < tolerance then (logLikelihood, parameters)
    else doStep parameters logLikelihood
  doStep initialParams -10e100
4

3 に答える 3

2

私は F# のやり方を推測しようとはしていませんが、いくつかの観察結果を以下に示します。

1) 観察した初期状態はいくつありますか? 答えが「1 つだけ」の場合、観測の確率は P(状態 0) P(obs | 状態は 0) + P(状態 1) P(obs | 状態 1) と書くことができます。2 つの P(obs | state is X) のどちらが高いかに応じて、最尤解は P(state 0) = 1 または P(state 1) = 1 になります。多数の異なる初期状態から派生する観測を観察している可能性がある場合の初期状態 - たとえば、同時に分析するおもちゃのデータの複数のストレッチがある場合。

2) バグを探す場合、答えが完全に明らかなおもちゃのデータを作成すると役立ちます。{0, 0, 0, 0...} という形式の n ストレッチのデータと、{1, 1, 1, 1...} という形式の m ストレッチのデータがある場合、状態 0 が割り当てられることを望むかもしれません。初期確率 n/(m +n) - またはもちろん m/(n + m) です。これは、プログラムがどの状態をどのシーケンスにリンクしたいかを知らないためです。

3) プログラムをチェックするもう 1 つの方法は、ある種の一貫性または保存チェックを探すことです。2 つの初期状態のモデルは、1 つの初期状態、最初の観測に対する特別な遷移確率のセット、および場合によっては特別なダミーの最初の観測を持つモデルと同じにすることができるため、2 つの初期状態での動作を次のようにチェックできます。初期状態が 1 つだけで、多少の操作が必要な場合の動作。

于 2012-08-08T17:37:34.027 に答える
0

しばらく考えてみたところ、私が説明した振る舞いは実際には正しいのではないかと思います。その理由は、観測データではこの分布が1回だけ「使用」されたため、直感的に分布を推測するのに十分な統計がないためです。とは言っても、アルゴリズムは、時系列全体に影響を与えたため、時間0での非表示状態変数の実際の値を（妥当な精度で）回復できるはずだと思います。

于 2012-08-09T09:15:16.087 に答える