本文主要是介绍C#，数值计算——隐式马尔科夫模型（Hidden Markov Models）的计算方法与源程序，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！
1 文本格式
using System;
namespace Legalsoft.Truffer
{
/// <summary>
/// Hidden Markov Models
/// </summary>
public class HMM
{
private int fbdone { get; set; }
private int mstat { get; set; }
private int nobs { get; set; }
private int ksym { get; set; }
private int lrnrm { get; set; }
private double BIG { get; set; }
private double BIGI { get; set; }
private double lhood { get; set; }
private double[,] a { get; set; }
private double[,] b { get; set; }
private int[] obs { get; set; }
private double[,] alpha { get; set; }
private double[,] beta { get; set; }
private double[,] pstate { get; set; }
private int[] arnrm { get; set; }
private int[] brnrm { get; set; }
public HMM(double[,] aa, double[,] bb, int[] obss)
{
this.a = aa;
this.b = bb;
this.obs = obss;
this.fbdone = 0;
this.mstat = a.GetLength(0);
this.nobs = obs.Length;
this.ksym = b.GetLength(1);
this.alpha = new double[nobs, mstat];
this.beta = new double[nobs, mstat];
this.pstate = new double[nobs, mstat];
this.arnrm = new int[nobs];
this.brnrm = new int[nobs];
this.BIG = 1.0e20;
this.BIGI = 1.0 / BIG;
if (a.GetLength(1) != mstat)
{
throw new Exception("transition matrix not square");
}
if (b.GetLength(0) != mstat)
{
throw new Exception("symbol prob matrix wrong size");
}
for (int i = 0; i < nobs; i++)
{
if (obs[i] < 0 || obs[i] >= ksym)
{
throw new Exception("bad data in obs");
}
}
for (int i = 0; i < mstat; i++)
{
double sum = 0.0;
for (int j = 0; j < mstat; j++)
{
sum += a[i, j];
}
if (Math.Abs(sum - 1.0) > 0.01)
{
throw new Exception("transition matrix not normalized");
}
for (int j = 0; j < mstat; j++)
{
a[i, j] /= sum;
}
}
for (int i = 0; i < mstat; i++)
{
double sum = 0.0;
for (int k = 0; k < ksym; k++)
{
sum += b[i, k];
}
if (Math.Abs(sum - 1.0) > 0.01)
{
throw new Exception("symbol prob matrix not normalized");
}
for (int k = 0; k < ksym; k++)
{
b[i, k] /= sum;
}
}
}
public double loglikelihood()
{
return Math.Log(lhood) + lrnrm * Math.Log(BIGI);
}
public void forwardbackward()
{
for (int i = 0; i < mstat; i++)
{
alpha[0, i] = b[i, obs[0]];
}
arnrm[0] = 0;
for (int t = 1; t < nobs; t++)
{
double asum = 0;
for (int j = 0; j < mstat; j++)
{
double sum = 0.0;
for (int i = 0; i < mstat; i++)
{
sum += alpha[t - 1, i] * a[i, j] * b[j, obs[t]];
}
alpha[t, j] = sum;
asum += sum;
}
arnrm[t] = arnrm[t - 1];
if (asum < BIGI)
{
++arnrm[t];
for (int j = 0; j < mstat; j++)
{
alpha[t, j] *= BIG;
}
}
}
for (int i = 0; i < mstat; i++)
{
beta[nobs - 1, i] = 1.0;
}
brnrm[nobs - 1] = 0;
for (int t = nobs - 2; t >= 0; t--)
{
double bsum = 0.0;
for (int i = 0; i < mstat; i++)
{
double sum = 0.0;
for (int j = 0; j < mstat; j++)
{
sum += a[i, j] * b[j, obs[t + 1]] * beta[t + 1, j];
}
beta[t, i] = sum;
bsum += sum;
}
brnrm[t] = brnrm[t + 1];
if (bsum < BIGI)
{
++brnrm[t];
for (int j = 0; j < mstat; j++)
{
beta[t, j] *= BIG;
}
}
}
lhood = 0.0;
for (int i = 0; i < mstat; i++)
{
lhood += alpha[0, i] * beta[0, i];
}
lrnrm = arnrm[0] + brnrm[0];
if (lhood != 0.0)
{
while (lhood < BIGI)
{
lhood *= BIG;
lrnrm++;
}
}
for (int t = 0; t < nobs; t++)
{
double sum = 0.0;
for (int i = 0; i < mstat; i++)
{
sum += (pstate[t, i] = alpha[t, i] * beta[t, i]);
}
// sum = lhood*pow(BIGI, lrnrm - arnrm[t] - brnrm[t]);
for (int i = 0; i < mstat; i++)
{
pstate[t, i] /= sum;
}
}
fbdone = 1;
}
public void baumwelch()
{
double[,] bnew = new double[mstat, ksym];
double[] powtab = new double[10];
for (int i = 0; i < 10; i++)
{
powtab[i] = Math.Pow(BIGI, i - 6);
}
if (fbdone != 1)
{
throw new Exception("must do forwardbackward first");
}
for (int i = 0; i < mstat; i++)
{
double denom = 0.0;
for (int k = 0; k < ksym; k++)
{
bnew[i, k] = 0.0;
}
for (int t = 0; t < nobs - 1; t++)
{
double term = (alpha[t, i] * beta[t, i] / lhood) * powtab[arnrm[t] + brnrm[t] - lrnrm + 6];
denom += term;
bnew[i, obs[t]] += term;
}
for (int j = 0; j < mstat; j++)
{
double num = 0.0;
for (int t = 0; t < nobs - 1; t++)
{
num += alpha[t, i] * b[j, obs[t + 1]] * beta[t + 1, j] * powtab[arnrm[t] + brnrm[t + 1] - lrnrm + 6] / lhood;
}
a[i, j] *= (num / denom);
}
for (int k = 0; k < ksym; k++)
{
bnew[i, k] /= denom;
}
}
b = bnew;
fbdone = 0;
}
/// <summary>
/// Markov Models and Hidden Markov Modeling
/// </summary>
/// <param name="atrans"></param>
/// <param name="xout"></param>
/// <param name="istart"></param>
/// <param name="seed"></param>
/// <exception cref="Exception"></exception>
public static void markovgen(double[,] atrans, int[] xout, int istart = 0, int seed = 1)
{
int m = atrans.GetLength(0);
int n = xout.Length;
//double[,] cum = new double[,](atrans);
double[,] cum = Globals.CopyFrom(atrans);
Ran ran = new Ran((ulong)seed);
if (m != atrans.GetLength(1))
{
throw new Exception("transition matrix must be square");
}
for (int i = 0; i < m; i++)
{
for (int ja = 1; ja < m; ja++)
{
cum[i, ja] += cum[i, ja - 1];
}
if (Math.Abs(cum[i, m - 1] - 1.0) > 0.01)
{
throw new Exception("transition matrix rows must sum to 1");
}
}
int j = istart;
xout[0] = j;
for (int ii = 1; ii < n; ii++)
{
double r = ran.doub() * cum[j, m - 1];
int ilo = 0;
int ihi = m;
while (ihi - ilo > 1)
{
int ia = (ihi + ilo) >> 1;
if (r > cum[j, ia - 1])
{
ilo = ia;
}
else
{
ihi = ia;
}
}
xout[ii] = j = ilo;
}
}
}
}
2 代码格式
using System;namespace Legalsoft.Truffer
{/// <summary>/// Hidden Markov Models/// </summary>public class HMM{private int fbdone { get; set; }private int mstat { get; set; }private int nobs { get; set; }private int ksym { get; set; }private int lrnrm { get; set; }private double BIG { get; set; }private double BIGI { get; set; }private double lhood { get; set; }private double[,] a { get; set; }private double[,] b { get; set; }private int[] obs { get; set; }private double[,] alpha { get; set; }private double[,] beta { get; set; }private double[,] pstate { get; set; }private int[] arnrm { get; set; }private int[] brnrm { get; set; }public HMM(double[,] aa, double[,] bb, int[] obss){this.a = aa;this.b = bb;this.obs = obss;this.fbdone = 0;this.mstat = a.GetLength(0);this.nobs = obs.Length;this.ksym = b.GetLength(1);this.alpha = new double[nobs, mstat];this.beta = new double[nobs, mstat];this.pstate = new double[nobs, mstat];this.arnrm = new int[nobs];this.brnrm = new int[nobs];this.BIG = 1.0e20;this.BIGI = 1.0 / BIG;if (a.GetLength(1) != mstat){throw new Exception("transition matrix not square");}if (b.GetLength(0) != mstat){throw new Exception("symbol prob matrix wrong size");}for (int i = 0; i < nobs; i++){if (obs[i] < 0 || obs[i] >= ksym){throw new Exception("bad data in obs");}}for (int i = 0; i < mstat; i++){double sum = 0.0;for (int j = 0; j < mstat; j++){sum += a[i, j];}if (Math.Abs(sum - 1.0) > 0.01){throw new Exception("transition matrix not normalized");}for (int j = 0; j < mstat; j++){a[i, j] /= sum;}}for (int i = 0; i < mstat; i++){double sum = 0.0;for (int k = 0; k < ksym; k++){sum += b[i, k];}if (Math.Abs(sum - 1.0) > 0.01){throw new Exception("symbol prob matrix not normalized");}for (int k = 0; k < ksym; k++){b[i, k] /= sum;}}}public double loglikelihood(){return Math.Log(lhood) + lrnrm * Math.Log(BIGI);}public void forwardbackward(){for (int i = 0; i < mstat; i++){alpha[0, i] = b[i, obs[0]];}arnrm[0] = 0;for (int t = 1; t < nobs; t++){double asum = 0;for (int j = 0; j < mstat; j++){double sum = 0.0;for (int i = 0; i < mstat; i++){sum += alpha[t - 1, i] * a[i, j] * b[j, obs[t]];}alpha[t, j] = sum;asum += sum;}arnrm[t] = arnrm[t - 1];if (asum < BIGI){++arnrm[t];for (int j = 0; j < mstat; j++){alpha[t, j] *= BIG;}}}for (int i = 0; i < mstat; i++){beta[nobs - 1, i] = 1.0;}brnrm[nobs - 1] = 0;for (int t = nobs - 2; t >= 0; t--){double bsum = 0.0;for (int i = 0; i < mstat; i++){double sum = 0.0;for (int j = 0; j < mstat; j++){sum += a[i, j] * b[j, obs[t + 1]] * beta[t + 1, j];}beta[t, i] = sum;bsum += sum;}brnrm[t] = brnrm[t + 1];if (bsum < BIGI){++brnrm[t];for (int j = 0; j < mstat; j++){beta[t, j] *= BIG;}}}lhood = 0.0;for (int i = 0; i < mstat; i++){lhood += alpha[0, i] * beta[0, i];}lrnrm = arnrm[0] + brnrm[0];if (lhood != 0.0){while (lhood < BIGI){lhood *= BIG;lrnrm++;}}for (int t = 0; t < nobs; t++){double sum = 0.0;for (int i = 0; i < mstat; i++){sum += (pstate[t, i] = alpha[t, i] * beta[t, i]);}// sum = lhood*pow(BIGI, lrnrm - arnrm[t] - brnrm[t]);for (int i = 0; i < mstat; i++){pstate[t, i] /= sum;}}fbdone = 1;}public void baumwelch(){double[,] bnew = new double[mstat, ksym];double[] powtab = new double[10];for (int i = 0; i < 10; i++){powtab[i] = Math.Pow(BIGI, i - 6);}if (fbdone != 1){throw new Exception("must do forwardbackward first");}for (int i = 0; i < mstat; i++){double denom = 0.0;for (int k = 0; k < ksym; k++){bnew[i, k] = 0.0;}for (int t = 0; t < nobs - 1; t++){double term = (alpha[t, i] * beta[t, i] / lhood) * powtab[arnrm[t] + brnrm[t] - lrnrm + 6];denom += term;bnew[i, obs[t]] += term;}for (int j = 0; j < mstat; j++){double num = 0.0;for (int t = 0; t < nobs - 1; t++){num += alpha[t, i] * b[j, obs[t + 1]] * beta[t + 1, j] * powtab[arnrm[t] + brnrm[t + 1] - lrnrm + 6] / lhood;}a[i, j] *= (num / denom);}for (int k = 0; k < ksym; k++){bnew[i, k] /= denom;}}b = bnew;fbdone = 0;}/// <summary>/// Markov Models and Hidden Markov Modeling/// </summary>/// <param name="atrans"></param>/// <param name="xout"></param>/// <param name="istart"></param>/// <param name="seed"></param>/// <exception cref="Exception"></exception>public static void markovgen(double[,] atrans, int[] xout, int istart = 0, int seed = 1){int m = atrans.GetLength(0);int n = xout.Length;//double[,] cum = new double[,](atrans);double[,] cum = Globals.CopyFrom(atrans);Ran ran = new Ran((ulong)seed);if (m != atrans.GetLength(1)){throw new Exception("transition matrix must be square");}for (int i = 0; i < m; i++){for (int ja = 1; ja < m; ja++){cum[i, ja] += cum[i, ja - 1];}if (Math.Abs(cum[i, m - 1] - 1.0) > 0.01){throw new Exception("transition matrix rows must sum to 1");}}int j = istart;xout[0] = j;for (int ii = 1; ii < n; ii++){double r = ran.doub() * cum[j, m - 1];int ilo = 0;int ihi = m;while (ihi - ilo > 1){int ia = (ihi + ilo) >> 1;if (r > cum[j, ia - 1]){ilo = ia;}else{ihi = ia;}}xout[ii] = j = ilo;}}}
}
这篇关于C#，数值计算——隐式马尔科夫模型（Hidden Markov Models）的计算方法与源程序的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！
