uhmm.pas

delphi 的hmm 程序

{@abstract(code for the tutorial on hidden Markov models )
@author(Nikolai Shokhirev <[e-mail]> <[e-mail]>
        http://www.shokhirev.com/nikolai.html
        http://www.u.arizona.edu/~nikolai/
        http://www.chem.arizona.edu/~shokhirn/nikolai.html )
@created(2006.02.02)
㎞ikolai V. Shokhirev, 2003-2006 }
{ 2006.02.24 - added PosteriorDecodingIdxs }
unit uHMM;

interface

uses
  uMatTypes, uDynArrays;

type

  Thmm = class(TObject)
  private
    N: TInt;
    M: TInt;
    Tmax: TInt;
    A: IFArr2D;
    B: IFArr2D;
    P0: IFArr1D;
    beta: IFArr2D;
    alpha: IFArr2D;
  public
//    procedure SetObservations(const ObsIdxs: IIArr1D);
    constructor Create(const NStates: TInt; const NObs: TInt); overload;
    constructor Create(const TransProb, ObsProb: IFArr2D;
                                     const InitProb: IFArr1D); overload;
    destructor Destroy; override;
    function GetProbabilityF(const ObsIdxs: IIArr1D): TFloat;
    function GetProbabilityB(const ObsIdxs: IIArr1D): TFloat;
    function ViterbiStateIdxs(const ObsIdxs: IIArr1D): IIArr1D;
    function PosteriorDecodingIdxs(const ObsIdxs: IIArr1D): IIArr1D;
    procedure SetTransitions(const TransProb: IFArr2D);
    procedure SetEmissions(const ObsProb: IFArr2D);
    procedure SetProbabilities(const InitProb: IFArr1D);
  end;

implementation

//uses


{ Thmm }

constructor Thmm.Create(const NStates, NObs: TInt);
begin
  inherited Create;
  N := NStates;
  M := NObs;
  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;
  A := TFArr2D.Create(TransProb, true);
  B := TFArr2D.Create(ObsProb, true);
  P0 := TFArr1D.Create(InitProb, true);
end;

destructor Thmm.Destroy;
begin
  A := nil;
  B := nil;
  P0 := nil;
  inherited;
end;

// Forward algorithm
function Thmm.GetProbabilityF(const ObsIdxs: IIArr1D): TFloat;
var
//  a1, a2: IFArr1D;
  t, i, j: TInt;
  sum: TFloat;
begin
//  a1 := TFArr1D.Create(A.Lo1, A.Hi1);
//  a2 := TFArr1D.Create(a1);
  alpha := TFArr2D.Create(ObsIdxs.Lo1, ObsIdxs.Hi1, A.Lo1, A.Hi1);

  // Initialization
  t := ObsIdxs.Lo1;
  for i := A.Lo1 to A.Hi1 do
    alpha[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 + alpha[t,i]*A[i,j];

      alpha[t+1,j] := sum*B[j,ObsIdxs[t+1]];
    end;
  end;

  // Termination
  t := ObsIdxs.Hi1;
  sum := 0.0;
  for i := A.Lo1 to A.Hi1 do
    sum := sum + alpha[t,i];

  result := sum;
end;

// Backward algorithm
function Thmm.GetProbabilityB(const ObsIdxs: IIArr1D): TFloat;
var
  t, i, j: TInt;
  sum: TFloat;
begin
  beta := TFArr2D.Create(ObsIdxs.Lo1, ObsIdxs.Hi1, A.Lo1, A.Hi1);

  // Initialization
  t := ObsIdxs.Hi1;
  for i := A.Lo1 to A.Hi1 do
    beta[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]]*beta[t+1,j];

      beta[t,i] := sum;
    end;
  end;

  // 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.ViterbiStateIdxs(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.PosteriorDecodingIdxs(const ObsIdxs: IIArr1D): IIArr1D;
var
  pB, pF: TFloat;
  t, i, k: TInt;
  p, q: TFloat;
  StateIdxs: IIArr1D;
begin
  StateIdxs := TIArr1D.Create(ObsIdxs.Lo1, ObsIdxs.Hi1);
  // calculate alpha
  pF := GetProbabilityF(ObsIdxs);
  // calculate beta
  pB := GetProbabilityB(ObsIdxs);


  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;

end.