Back
{@abstract(code for the tutorial on hidden Markov models http://www.shokhirev.com/nikolai/abc/alg/hmm/hmm.html )
@author(Nikolai Shokhirev <nikolai@shokhirev.com> http://www.shokhirev.com/nikolai.html
@created(2006.02.02)
ŠNikolai V. Shokhirev, 2003-2006 }
{ 2006.02.24 - added PosteriorDecodingIdxs }
{ 2009.05.16 - refactored ForwardProbability, BackwardProbability; added GetAlpha, GetBeta }
{ 2009.05.17 - added BaumWelchIteration, GetNorm, GetGamma, GetP0Estimate, GetAEstimate, GetBEstimate}
{ 2010.02.15 - Cleanup}
unit uHMM;
interface
uses
uMatTypes, uDynArrays;
type
Thmm = class(TObject)
private
N: TInt; // Number of states
M: TInt; // Number of observation
A: IFArr2D; // Transition probabilities A[1..N, 1..N]
B: IFArr2D; // Observation probabilities B[1..N, 1..M]
P0: IFArr1D; // Initial probabilities P0[1..N]
// Tmax - Number of observations
// ObsIdxs: IIArr1D - Indices of observations ObsIdxs[1..Tmax]
// returns Forward probability array [1..Tmax, 1..N]
function GetAlpha(const ObsIdxs: IIArr1D; const A, B: IFArr2D;
const P0: IFArr1D): IFArr2D;
// returns Backward probability array [1..Tmax, 1..N]
function GetBeta(const ObsIdxs: IIArr1D; const A, B: IFArr2D): IFArr2D;
// returns Sum(i = 1 to N, alpha[t,i]*beta[t,i]) - Does not dpend on t
function GetNorm(const alpha, beta: IFArr2D): TFloat;
// returns gamma [1..Tmax, 1..N]
function GetGamma(const alpha, beta: IFArr2D; Norm: TFloat): IFArr2D;
function GetP0Estimate(const gamma: IFArr2D): IFArr1D;
function GetAEstimate(const ObsIdxs: IIArr1D; const A, B: IFArr2D;
const alpha, beta, gamma: IFArr2D;
Norm: TFloat): IFArr2D;
function GetBEstimate(const ObsIdxs: IIArr1D; const gamma: IFArr2D): IFArr2D;
public
constructor Create(const NStates: TInt; const NObs: TInt); overload;
constructor Create(const TransProb, ObsProb: IFArr2D;
const InitProb: IFArr1D); overload;
destructor Destroy; override;
function ForwardProbability(const ObsIdxs: IIArr1D): TFloat;
function BackwardProbability(const ObsIdxs: IIArr1D): TFloat;
function ViterbiIndices(const ObsIdxs: IIArr1D): IIArr1D;
function PosteriorIndices(const ObsIdxs: IIArr1D): IIArr1D;
procedure SetTransitions(const TransProb: IFArr2D);
procedure SetEmissions(const ObsProb: IFArr2D);
procedure SetProbabilities(const InitProb: IFArr1D);
procedure BaumWelchIteration(const ObsIdxs: IIArr1D;
var NewA, NewB: IFArr2D; var NewP0: IFArr1D);
end;
implementation
uses SysUtils, Dialogs;
{ Thmm }
constructor Thmm.Create(const NStates, NObs: TInt);
begin
inherited Create;
N := NStates;
M := NObs;
// reference-counting arrays; automatically destroyed when go out of scope.
// Defined in uDynArrays. Download http://www.shokhirev.com/nikolai/programs/tools/PasMatLib/download.html .
A := TFArr2D.Create(N,N);
B := TFArr2D.Create(N,M);
P0 := TFArr1D.Create(N);
end;
constructor Thmm.Create(const TransProb, ObsProb: IFArr2D;
const InitProb: IFArr1D);
begin
inherited Create;
if (not SameLim1(TransProb, ObsProb)) or (not SameLim1(TransProb, InitProb)) then
raise ELimMismatch.Create('Lim1 mismatch for TransProb, ObsProb, InitProb');
if not IsSquare(TransProb) then
raise ENonSquareMatrix.Create('TransProb is not square');
N := TransProb.Dim1;
M := ObsProb.Dim2;
// reference-counting arrays; automatically destroyed when go out of scope.
A := TFArr2D.Create(TransProb, true); // Clone
B := TFArr2D.Create(ObsProb, true); // Clone
P0 := TFArr1D.Create(InitProb, true); // Clone
end;
destructor Thmm.Destroy;
begin
A := nil;
B := nil;
P0 := nil;
inherited;
end;
// returns Forward probability array [1..Tmax, 1..N]
function Thmm.GetAlpha(const ObsIdxs: IIArr1D; const A: IFArr2D;
const B: IFArr2D; const P0: IFArr1D): IFArr2D;
var
t, i, j: TInt;
sum: TFloat;
begin
result := TFArr2D.Create(ObsIdxs.Lo1, ObsIdxs.Hi1, A.Lo1, A.Hi1);
// Initialization
t := ObsIdxs.Lo1;
for i := A.Lo1 to A.Hi1 do
result[t,i] := P0[i]*B[i,ObsIdxs[t]];
// Recursion
for t := ObsIdxs.Lo1 to ObsIdxs.Hi1-1 do
begin
for j := A.Lo1 to A.Hi1 do
begin
sum := 0.0;
for i := A.Lo1 to A.Hi1 do
sum := sum + result[t,i]*A[i,j];
result[t+1,j] := sum*B[j,ObsIdxs[t+1]];
end;
end;
end;
// Forward algorithm
function Thmm.ForwardProbability(const ObsIdxs: IIArr1D): TFloat;
var
t, i: TInt;
sum: TFloat;
alpha: IFArr2D; // Forward probabilities alpha.Dim1 = Tmax, alpha.Dim2 = N
begin
// Initialization, Recursion
alpha := GetAlpha(ObsIdxs, A, B, P0);
// Termination
t := ObsIdxs.Hi1;
sum := 0.0;
for i := A.Lo1 to A.Hi1 do
sum := sum + alpha[t,i];
result := sum;
end;
// returns Backward probability array [1..Tmax, 1..N]
function Thmm.GetBeta(const ObsIdxs: IIArr1D; const A: IFArr2D;
const B: IFArr2D): IFArr2D;
var
t, i, j: TInt;
sum: TFloat;
begin
result := TFArr2D.Create(ObsIdxs.Lo1, ObsIdxs.Hi1, A.Lo1, A.Hi1);
// Initialization
t := ObsIdxs.Hi1;
for i := A.Lo1 to A.Hi1 do
result[t,i] := 1.0;
// Recursion
for t := ObsIdxs.Hi1-1 downto ObsIdxs.Lo1 do
begin
for i := A.Lo1 to A.Hi1 do
begin
sum := 0.0;
for j := A.Lo1 to A.Hi1 do
sum := sum + A[i,j]*B[j,ObsIdxs[t+1]]*result[t+1,j];
result[t,i] := sum;
end;
end;
end;
// Backward algorithm
function Thmm.BackwardProbability(const ObsIdxs: IIArr1D): TFloat;
var
t, i: TInt;
sum: TFloat;
beta: IFArr2D; // Backward probabilities beta.Dim1 = Tmax, beta.Dim2 = N
begin
// Initialization, Recursion
beta := GetBeta(ObsIdxs, A, B);
// Termination
t := ObsIdxs.Lo1;
sum := 0.0;
for i := A.Lo1 to A.Hi1 do
sum := sum + P0[i]*B[i,ObsIdxs[t]]*beta[t,i];
result := sum;
end;
// Viterbi algorithm
function Thmm.ViterbiIndices(const ObsIdxs: IIArr1D): IIArr1D;
var
StateIdxs: IIArr1D;
delta: IFArr2D;
psi: IIArr2D;
t, i, j, k: TInt;
p, q: TFloat;
begin
StateIdxs := TIArr1D.Create(ObsIdxs.Lo1, ObsIdxs.Hi1);
delta := TFArr2D.Create(ObsIdxs.Lo1, ObsIdxs.Hi1, A.Lo1, A.Hi1);
psi := TIArr2D.Create(ObsIdxs.Lo1, ObsIdxs.Hi1, A.Lo1, A.Hi1);
// Initialization
t := ObsIdxs.Lo1;
for i := A.Lo1 to A.Hi1 do
begin
delta[t,i] := P0[i]*B[i,ObsIdxs[t]];
psi[t,i] := A.Lo1 - 1; // Outside limits - not used
end;
// Recursion
for t := ObsIdxs.Lo1+1 to ObsIdxs.Hi1 do
begin
for j := A.Lo1 to A.Hi1 do
begin
p := 0.0;
k := A.Lo1-1; // Outside limits, must be redefined below
for i := A.Lo1 to A.Hi1 do
begin
q := delta[t-1,i]*A[i,j];
if q > p then
begin
p := q;
k := i;
end;
end;
delta[t,j] := p*B[j,ObsIdxs[t]];
psi[t,j] := k;
end;
end;
// Termination
t := ObsIdxs.Hi1;
p := 0.0;
k := A.Lo1-1; // Outside limits, must be redefined below
for i := A.Lo1 to A.Hi1 do
begin
q := delta[t,i];
if q > p then
begin
p := q;
k := i;
end;
end;
StateIdxs[t] := k; // q* in Rabiner's paper
// Path (state sequence) backtracking
for t := ObsIdxs.Hi1-1 downto ObsIdxs.Lo1 do
begin
StateIdxs[t] := psi[t+1, StateIdxs[t+1]];
end;
result := StateIdxs;
end;
// Posterior Decoding
function Thmm.PosteriorIndices(const ObsIdxs: IIArr1D): IIArr1D;
var
// pB, pF: TFloat;
t, i, k: TInt;
p, q: TFloat;
StateIdxs: IIArr1D;
alpha: IFArr2D; // Forward probabilities alpha.Dim1 = Tmax, alpha.Dim2 = N
beta: IFArr2D; // Backward probabilities beta.Dim1 = Tmax, beta.Dim2 = N
begin
StateIdxs := TIArr1D.Create(ObsIdxs.Lo1, ObsIdxs.Hi1);
// calculate alpha
alpha := GetAlpha(ObsIdxs, A, B, P0);
// calculate beta
beta := GetBeta(ObsIdxs, A, B);
for t := ObsIdxs.Lo1 to ObsIdxs.Hi1 do
begin
p := 0.0;
k := A.Lo1-1; // Outside limits, must be redefined below
for i := A.Lo1 to A.Hi1 do
begin
q := alpha[t,i]*beta[t,i];
if q > p then
begin
p := q;
k := i;
end;
end;
StateIdxs[t] := k;
end;
result := StateIdxs;
end;
procedure Thmm.SetEmissions(const ObsProb: IFArr2D);
begin
if (not SameLim(B, ObsProb)) then
raise ELimMismatch.Create('Lim mismatch for B, ObsProb');
B.Assign(ObsProb);
end;
procedure Thmm.SetProbabilities(const InitProb: IFArr1D);
begin
if (not SameLim(P0, InitProb)) then
raise ELimMismatch.Create('Lim1 mismatch for P0, InitProb');
P0.Assign(InitProb);
end;
procedure Thmm.SetTransitions(const TransProb: IFArr2D);
begin
if (not SameLim(A, TransProb)) then
raise ELimMismatch.Create('Lim mismatch for A, TransProb');
A.Assign(TransProb);
end;
// result = Sum(i = 1 to N, alpha[t,i]*beta[t,i]) - Does not dpend on t
// = ForwardProbability = BackwardProbability
function Thmm.GetNorm(const alpha, beta: IFArr2D): TFloat;
var
t, i: TInt;
sum: TFloat;
begin
t := alpha.Lo1;
sum := 0.0;
for i := alpha.Lo2 to alpha.Hi2 do
sum := sum + alpha[t,i]*beta[t,i];
result := sum;
end;
// returns gamma [1..Tmax, 1..N]
function Thmm.GetGamma(const alpha, beta: IFArr2D; Norm: TFloat): IFArr2D;
var
t, i: TInt;
begin
result := TFArr2D.Create(alpha);
for t := alpha.Lo1 to alpha.Hi1 do
begin
for i := alpha.Lo2 to alpha.Hi2 do
result[t,i] := alpha[t,i]*beta[t,i]/Norm;
end;
end;
function Thmm.GetP0Estimate(const gamma: IFArr2D): IFArr1D;
var
t, i: TInt;
begin
result := TFArr1D.Create(gamma.Lo2, gamma.Hi2);
t := gamma.Lo1;
for i := gamma.Lo2 to gamma.Hi2 do
result[i] := gamma[t,i];
end;
function Thmm.GetAEstimate(const ObsIdxs: IIArr1D; const A, B: IFArr2D;
const alpha, beta, gamma: IFArr2D;
Norm: TFloat): IFArr2D;
var
t, i, j: TInt;
num, den, xi: TFloat;
begin
result := TFArr2D.Create(A);
for i := A.Lo1 to A.Hi1 do
begin
den := 0.0;
for t := ObsIdxs.Lo1 to ObsIdxs.Hi1-1 do
den := den + gamma[t,i];
for j := A.Lo2 to A.Hi2 do
begin
num := 0.0;
for t := ObsIdxs.Lo1 to ObsIdxs.Hi1-1 do
begin
xi := alpha[t,i]*A[i,j]*B[j,ObsIdxs[t+1]]*beta[t+1,i]/Norm;
num := num + xi;
end;
result[i,j] := num/den;
end;
end;
end;
function Thmm.GetBEstimate(const ObsIdxs: IIArr1D; const gamma: IFArr2D): IFArr2D;
var
t, i, k: TInt;
den: TFloat;
begin
result := TFArr2D.Create(B);
result.Fill(0.0);
for i := A.Lo1 to A.Hi1 do
begin
den := 0.0;
for t := ObsIdxs.Lo1 to ObsIdxs.Hi1 do
den := den + gamma[t,i];
for t := ObsIdxs.Lo1 to ObsIdxs.Hi1 do
begin
k := ObsIdxs[t];
result[i,k] := result[i,k] + gamma[t,i]/den;
end;
end;
end;
procedure Thmm.BaumWelchIteration(const ObsIdxs: IIArr1D;
var NewA, NewB: IFArr2D; var NewP0: IFArr1D);
var
alpha: IFArr2D; // Forward probabilities alpha.Dim1 = Tmax, alpha.Dim2 = N
beta: IFArr2D; // Backward probabilities beta.Dim1 = Tmax, beta.Dim2 = N
gamma: IFArr2D; // gamma.Dim1 = Tmax, gamma.Dim2 = N
Norm: TFloat;
begin
alpha := GetAlpha(ObsIdxs, A, B, P0);
beta := GetBeta(ObsIdxs, A, B);
// returns Sum(i = 1 to N, alpha[t,i]*beta[t,i]) - Does not dpend on t
Norm := GetNorm(alpha, beta);
gamma := GetGamma(alpha, beta, Norm);
// estimates
NewP0 := GetP0Estimate(gamma);
NewA := GetAEstimate(ObsIdxs, A, B, alpha, beta, gamma, Norm);
NewB := GetBEstimate(ObsIdxs, gamma, Norm);
end;
end.