HMM model 例子 Biased coins

最新推荐文章于 2024-08-08 15:12:59 发布
最新推荐文章于 2024-08-08 15:12:59 发布
阅读量1.2k 收藏
点赞数
分类专栏： EM HMM 文章标签： EM 算法 HMM C++
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/u012551177/article/details/18638809
版权


输入

$ head toss1.txt
H
T
H
T
T
H
T
H
T
T
H for 正面， T for 反面。
可能是正常的 coin，可能是biased coin（HT 的概率非0.5 0.5），二者按照一定几率转换，几率不知道。


输出
TIME    TOSS    P(FAIR) P(BIAS) MLSTATE
1       H       0.3131  0.6869  BIASED
2       T       0.9991  0.0009  FAIR
3       H       0.9788  0.0212  FAIR
4       T       0.9978  0.0022  FAIR
5       T       0.9978  0.0022  FAIR
6       H       0.9788  0.0212  FAIR
7       T       0.9978  0.0022  FAIR
8       H       0.9788  0.0212  FAIR
9       T       0.9978  0.0022  FAIR
10      T       0.9978  0.0022  FAIR


根据 sph.umich  Biostatistics 615课程课件改编。


main.cpp

这里假定训练参数到 0.001截止，否则会过度训练，最多通过 EM算法迭代100次。

#include <iostream>
#include <iomanip>
#include "HMM615.h"

using namespace std;

int main(int argc, char** argv) {
	double tran_threshold = 0.001;
	vector<int> toss; string tok;
	while( cin >> tok ) {
		if ( tok == "H" ) toss.push_back(0);
		else if ( tok == "T" ) toss.push_back(1);
		else {
			cerr << "Cannot recognize input " << tok << endl;
			return -1;
		}
	}

	int T = toss.size();
	HMM615 hmm(2, 2, T);
	hmm.trans.data[0][0] = 0.95; hmm.trans.data[0][1] = 0.05;
	hmm.trans.data[1][0] = 0.2; hmm.trans.data[1][1] = 0.8;
	hmm.emis.data[0][0] = 0.5; hmm.emis.data[0][1] = 0.5;
	hmm.emis.data[1][0] = 0.9; hmm.emis.data[1][1] = 0.1;
	hmm.pis[0] = 0.5; hmm.pis[1] = 0.5; hmm.outs = toss;
	hmm.viterbi();
	hmm.BaumWelch();

	for (int i=1; i<100; ++i){
		if ( hmm.TRANS.data[0][0] < tran_threshold || hmm.TRANS.data[0][1] < tran_threshold || hmm.TRANS.data[1][0] < tran_threshold || hmm.TRANS.data[1][1] < tran_threshold || hmm.EMIS.data[0][0] < tran_threshold || hmm.EMIS.data[0][1] < tran_threshold || hmm.EMIS.data[1][0] < tran_threshold || hmm.EMIS.data[1][1] < tran_threshold  ){
			break;
		}
		hmm.trans.data[0][0] = hmm.TRANS.data[0][0];
  		hmm.trans.data[0][1] = hmm.TRANS.data[0][1];
		hmm.trans.data[1][0] = hmm.TRANS.data[1][0];
		hmm.trans.data[1][1] = hmm.TRANS.data[1][1];

		hmm.emis.data[0][0] = hmm.EMIS.data[0][0];
		hmm.emis.data[0][1] = hmm.EMIS.data[0][1];
		hmm.emis.data[1][0] = hmm.EMIS.data[1][0];
		hmm.emis.data[1][1] = hmm.EMIS.data[1][1];

		hmm.zero_sigma();
		hmm.viterbi();
		hmm.BaumWelch();
	}
	cout << "TIME\tTOSS\tP(FAIR)\tP(BIAS)\tMLSTATE" << endl;
	cout << setiosflags(ios::fixed) << setprecision(4);
	for(int t=0; t < T; ++t) {
		cout << t+1 << "\t" << (toss[t] == 0 ? "H" : "T") << "\t" << hmm.gammas.data[t][0] << "\t" << hmm.gammas.data[t][1] << "\t" << (hmm.path[t] == 0 ? "FAIR" : "BIASED" ) << endl;
	}
	cout << hmm.TRANS.data[0][0] << "\t" << hmm.TRANS.data[0][1] << endl;
	cout << hmm.TRANS.data[1][0] << "\t" << hmm.TRANS.data[1][1] << endl;
	cout << hmm.EMIS.data[0][0]  << "\t" << hmm.EMIS.data[0][1]  << endl;
	cout << hmm.EMIS.data[1][0]  << "\t" << hmm.EMIS.data[1][1]  << endl;
	return 0;
}

HMM615.h  

#ifndef __HMM_615_H
#define __HMM_615_H
#include "Matrix615.h"
#include "MatrixTrible.h"
#include <cmath>

class HMM615 {
	public:
	// parameters
		int nStates; // n : number of possible states
		int nObs; // m : number of possible output values
		int nTimes; // t : number of time slots with observations
		std::vector<double> pis; // initial states
		std::vector<int> outs; // observed outcomes
		Matrix615<double> trans; // trans[i][j] corresponds to A_{ij}
		Matrix615<double> emis;
		Matrix615<double> TRANS;	// Revised trans and emis
		Matrix615<double> EMIS;	//

	// storages for dynamic programming
	Matrix615<double> alphas, betas, gammas, deltas;
	Matrix615<int> phis;
	MatrixTrible<double> sigma;
	std::vector<int> path;
	HMM615(int states, int obs, int times) : nStates(states), nObs(obs), nTimes(times), trans(states, states, 0), emis(states, obs, 0), alphas(times, states, 0), betas(times, states, 0),gammas(times, states, 0), deltas(times, states, 0),sigma(times,states,states,0),phis(times, states, 0),TRANS(times,states,0),EMIS(times,states,0)
	{
	   pis.resize(nStates);
	   path.resize(nTimes);
	}
	void forward(); // given below
	void backward(); //
	void forwardBackward(); // given below
	void BaumWelch();
	void zero_sigma();
	void viterbi(); //
};
#endif // __HMM_615_H

void HMM615::zero_sigma(){
	for(int t=0; t < nTimes-1; ++t) {
		for(int i=0; i < nStates; ++i){
			for (int j=0; j < nStates; ++j){
				sigma.data[t][i][j] = 0;
			}
		}
	}
}

void HMM615::forward(){
   for(int i=0; i < nStates; ++i){
		double tmp = std::log(pis[i]) + std::log(emis.data[i][outs[0]]);
      alphas.data[0][i] = std::exp(tmp);
   }
   for(int t=1; t < nTimes; ++t){
      for(int i=0; i < nStates; ++i){
         alphas.data[t][i] = 0;
         for(int j=0; j < nStates; ++j) {
				double tmp = std::log( alphas.data[t-1][j] ) + std::log( trans.data[j][i] ) + std::log( emis.data[i][outs[t]]) ;
            alphas.data[t][i] += std::exp(tmp);
         }
      }
   }
}

void HMM615::backward() {
	for(int i=0; i < nStates; ++i) {
		betas.data[nTimes-1][i] = 1;
	}
	for(int t=nTimes-2; t >=0; --t) {
		for(int i=0; i < nStates; ++i) {
			betas.data[t][i] = 0;
			for(int j=0; j < nStates; ++j) {
				double tmp = std::log( betas.data[t+1][j] ) + std::log( trans.data[i][j] ) + std::log( emis.data[j][outs[t+1]]) ;
				betas.data[t][i] += std::exp(tmp);
			}
		}
	}
}

void HMM615::forwardBackward() {
	forward();
	backward();
	for(int t=0; t < nTimes; ++t) {
		double sum = 0;
		for(int i=0; i < nStates; ++i) {
			double tmp = std::log( alphas.data[t][i] ) + std::log( betas.data[t][i] );
			sum += std::exp(tmp);
		}
		for(int i=0; i < nStates; ++i) {
			double tmp = std::log( alphas.data[t][i] ) + std::log( betas.data[t][i] ) - std::log( sum );
			gammas.data[t][i] = std::exp(tmp);
		}
	}
	for(int t=0; t < nTimes-1; ++t) {
		double sum = 0;
		for(int i=0; i < nStates; ++i){
			for (int j=0; j < nStates; ++j){
				double tmp = std::log( alphas.data[t][i] ) + std::log( trans.data[i][j] ) + std::log( emis.data[j][outs[t+1]] ) + std::log( betas.data[t+1][j] );
				sum += std::exp(tmp);
			}
		}
		for(int i=0; i < nStates; ++i){
			for (int j=0; j < nStates; ++j){
				double tmp = std::log( alphas.data[t][i] ) + std::log( trans.data[i][j] ) + std::log( emis.data[j][outs[t+1]] ) + std::log( betas.data[t+1][j] ) - std::log( sum );
				sigma.data[t][i][j] += std::exp(tmp);
			}
		}
	}
}

void HMM615::BaumWelch() {
	forwardBackward();
	for (int i=0; i<nStates; ++i){
		double sum_gamma1 = 0;
		for (int t=0; t<nTimes-1; ++t){
			sum_gamma1 += gammas.data[t][i];
		}
		for (int j=0; j<nStates; ++j){
			double sum_sigma1 = 0;
			for (int t=0; t<nTimes-1; ++t){
				sum_sigma1 += sigma.data[t][i][j];
			}
			TRANS.data[i][j] = sum_sigma1 / sum_gamma1;
		}
	}

	for (int i=0; i<nStates; ++i){
		for (int k=0; k<nObs; ++k){
			double sum_gamma2 = 0;
			double sum_gamma3 = 0;
			for (int t=0; t<nTimes; ++t){
				sum_gamma2 += gammas.data[t][i];
				if (outs[t]==k){
					sum_gamma3 += gammas.data[t][i];
				}
			}
			EMIS.data[i][k]= sum_gamma3 / sum_gamma2;
		}
	}
}

void HMM615::viterbi() {
	for(int i=0; i < nStates; ++i) {
		deltas.data[0][i] = pis[i] * emis.data[i][ outs[0] ];
	}
	for(int t=1; t < nTimes; ++t) {
		for(int i=0; i < nStates; ++i) {
			int maxIdx = 0;
			double tmp = std::log( deltas.data[t-1][0] ) + std::log( trans.data[0][i] ) + std::log( emis.data[i][ outs[t] ] );
			double maxVal = std::exp(tmp);
			for(int j=1; j < nStates; ++j) {
				double tmp2 = std::log( deltas.data[t-1][j] ) + std::log( trans.data[j][i] ) + std::log( emis.data[i][ outs[t] ] );
				double val  = std::exp( tmp2 );
				if ( val > maxVal ) {
					maxIdx = j;
					maxVal = val;
				}
			}
			deltas.data[t][i] = maxVal;
			phis.data[t][i] = maxIdx;
		}
	}
	double maxDelta = deltas.data[nTimes-1][0];
	path[nTimes-1] = 0;
	for(int i=1; i < nStates; ++i){
		if ( maxDelta < deltas.data[nTimes-1][i] ) {
			maxDelta = deltas.data[nTimes-1][i];
			path[nTimes-i] = i;
		}
	}
	for(int t=nTimes-2; t >= 0; --t){
		path[t] = phis.data[t+1][ path[t+1] ];
	}
}

Matrix615.h 

#ifndef __MATRIX_615_H
#define __MATRIX_615_H
#include <vector>

template <class T>
class Matrix615 {
   public:
      std::vector< std::vector<T> > data;
      Matrix615(int nrow, int ncol, T val = 0) {
         data.resize(nrow); // make n rows
         for(int i=0; i < nrow; ++i) {
            data[i].resize(ncol,val); // make n cols with default value val
         }
      }
      int rowNums() { return (int)data.size(); }
      int colNums() { return ( data.size() == 0 ) ? 0 : (int)data[0].size(); }
};
#endif // __MATRIX_615_H

MatrixTrible.h 

#ifndef __MATRIX_TRIBLE_H
#define __MATRIX_TRIBLE_H
#include <vector>

template <class T>
class MatrixTrible {
   public:
      std::vector< std::vector< std::vector<T> > > data;
      MatrixTrible(int nrow, int ncol, int ntri, T val = 0) {
         data.resize(nrow); // make n rows
         for(int i=0; i < nrow; ++i) {
            data[i].resize(ncol); // make n cols with default value val
				for (int j=0; j < ncol; ++j){
					data[i][j].resize(ntri,val);
				}
         }
      }
      int rowNums() { return (int)data.size(); }
      int colNums() { return ( data.size() == 0 ) ? 0 : (int)data[0].size(); }
      int triNums() { return ( data.size() == 0 ) ? 0 : ( ( data[0].size() == 0) ? 0 : (int)data[0][0].size() ); }
};
#endif // __MATRIX_TRIBLE_H
Makefile 
cc=g++
obj=main.o
exe=hmm
$(exe):$(obj)
	$(cc) -o $(exe) $(obj)

main.o:main.cpp HMM615.h
	$(cc) -c main.cpp

clean:
	rm -rf *.o $(exe)



biopic
关注
  • 0
    点赞
  • 踩
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
通俗理解隐马尔可夫模型（HMM）及其案例代码实现
David's Tweet
05-08 4580
必要的数学知识 1.联合概率与边缘概率 联合概率是指多维随机变量中同时满足多个变量时候的概率，也就是共同发生的概率。A,B的联合概率通常写成P(A∩B)或P(AB)或P(AB)。对于离散的变量，联合概率可以用表格形式表示或者求和表示，连续的变量可以使用积分表示（若是二维就一个二重积分） 边缘概率是指多维随机变量中只满足部分变量时的概率。图片帮助理解： 联合概率与边缘概率的关系...
python翻转_Python程序设计偏向硬币翻转功能
cumubi7552的博客
07-04 531
python翻转Here, we are going to build a biasedcoin() function using python. The program is so simple as an introductory program and similar to the function coin() for defining a biased coin flip. The fu...
抛硬币正面期望_如果抛硬币，正面的数量多于反面的可能性
cumtb2002的博客
06-23 1457
抛硬币正面期望Problem statement: 问题陈述： Let N be a positive odd number. There are N coins, numbered 1, 2 , ... , N. For each i (1≤i≤N), when Coin i is tossed, it comes up heads with probability pi and tails...
斯坦福 机器学习-第四章学习理论
hzp666的博客
08-23 1601
CS229 Lecture notes 原作者：Andrew Ng（吴恩达） 翻译：CycleUser Part V
学习理论（Learning Theory） 1 偏差/方差的权衡（Bias/variance tradeoff ） 在讲线性回归的时候，我们讨论过这样的问题：拟合数据的时候，选择线性的“y = θ0 +θ1x”这样的“简单”模型，还是选择多项式
随机概率算法
maxam0128的专栏
01-07 1276
转载：http://www.keithschwarz.com/interesting/code/?dir=alias-method /******************************************************************************  * File: AliasMethod.java  * Author: Keith Schwarz 
HMM模型
zakexu的专栏
02-19 3215
（一）基本概念 1.隐马尔科夫模型（Hidden Markov Model，HMM）：关于时序的概率模型；描述由一个隐藏的马尔科夫链随机生成不可观测的状态随机序列，再由各个状态生成一个观测而产生观测随机序列的过程； （1）状态序列：HMM生成的状态的序列，称为状态序列； （2）观测序列：每个状态生成一个观测，而由此产生的观测的随机序列，称为观测序列； （3）序列的每一个位置又可以看作是一个...
HMM隐马尔科夫模型与实例1
sanjianjixiang的博客
12-26 639
隐马尔可夫模型 Hidden Markov Model (HMM) 以下三个问题，人们提出了相应的算法 1 评估问题： 前向算法 2 解码问题： Viterbi算法 (维特比算法) 3 学习问题： Baum-Welch算法(向前向后算法) import numpy as np from hmmlearn import hmm # 隐藏状态：3个盒子 states = ['box1', 'bo...
hmm.zip_HMM_generative model_vectors
09-21
在给定的"**hmm.zip_HMM_generative_model_vectors**"中，我们可以看到一个名为"hmm.m"的文件，这很可能是用MATLAB编写的一个实现HMM生成模型的脚本或函数。 隐藏马尔可夫模型（HMM）的核心概念是“隐藏状态”和...
HMM_model.zip_C HMM_HMM
09-22
隐马尔可夫模型（Hidden Markov Model，简称HMM）是统计建模方法中的一种，广泛应用于自然语言处理、语音识别、生物信息学等领域。HMM的核心思想是假设观察序列是由不可见的状态序列生成的，而状态之间的转移以及...
HMM.rar_HMM_MARKOV_in_markov model csharp
09-23
隐藏马尔可夫模型（Hidden Markov Model，简称HMM）是概率统计模型，它在计算机科学、机器学习和自然语言处理等领域有着广泛的应用。在C#编程语言中实现HMM，可以为各种任务提供强大的工具，比如语音识别、文本分类...
HMM.rar_ hmm_HMM例子_hmm mfc_hmm source codes_识别
09-22
隐马尔可夫模型（Hidden Markov Model，简称HMM）是统计建模方法中的一种，广泛应用于自然语言处理、语音识别、生物信息学等领域。HMM的核心思想是假设观察序列是由不可见的状态序列生成的，而状态之间的转移以及...
HMM.rar_HMM_hidden markov model _隐马尔科夫_马尔科夫
09-22
**隐马尔科夫模型（Hidden Markov Model，HMM）**是一种统计建模方法，广泛应用于自然语言处理、语音识别、生物信息学等领域。它基于马尔科夫假设，即系统状态的转移只依赖于当前状态，而与之前的状态无关。在HMM中...
CSP 2020 第三题：表达式
m0_46183615的博客
08-03 545
CSP 2020 第三题：表达式
VRF 算法流程及代码实现
smilejiasmile的博客
08-03 233
Verifiable Random Functions (VRFs) 是一种密码学原语，能够为给定输入生成一个随机输出，并提供一个可以验证的证明，证明输出是由特定输入生成的。VRFs 具有许多应用，包括区块链、分布式系统和安全协议。使用secp256k1曲线可以产生确定性签名，适用于 VRF 实现。通过使用 HMAC-SHA256 生成确定性随机数k，可以确保签名过程对相同的输入始终产生相同的结果，从而满足 VRF 的要求。
day 9 回溯算法
最新发布
m0_73924943的博客
08-08 145
优化后的比之前的在叶子结点少递归一层。在排序candidates后，碰到target<0的情况直接结束本层的遍历，后面没遍历到的也不用遍历了。1.目标数n一直在随i的加入减小，如果目标数n已经<0，那么剪枝。即path中的数已经大于目标数了。利用数组去重，因为题目的数值范围小，所以可以尽量使用。剪枝优化：（对每层循环遍历的范围进行了缩小）2.对for循环i的遍历范围进行缩小。也可以不用st数组去重。
嵌入式初学-C语言-十五
MaBaB的博客
08-05 616
函数：实现一定功能的，独立的代码模块，我们的函数一定是先定义后使用使用函数的优势可以通过函数提供功能给别人使用，当然我们也可以使用别人提供的函数，减少代码量借助函数可以减少重复性代码。实现结构化（模块化）程序设计思想结构化程序设计思想：将大型的任务功能，划分为相互独立的小型的任务模块来设计3.函数是c语言程序的基本组成单元：C语言程序是由一个（必然是main函数）或多个函数组成。语法：返回类型 函数名(形参列表) -- 函数头 | 函数首部。
【算法】动态规划解决背包问题
2302_78914800的博客
08-04 456
1.动态规划算法的核心是：将大问题划分为小问题进行解决，从而一步步获取最优解的处理算法2.动态规划算法与分治算法类似，其基本思想也是将待求解问题分解成若干个子问题，先求解子问题，然后从这些子问题的解得到原问题的解3.与分治算法不同的是，适合于用动态规划求解的问题，经分解得到的子问题往往不是相互独立的（即下一个子阶段的求解是建立在上一个子阶段的解的基础上，进行进一步的求解）4.动态规划可以通过填表的方式来逐步推进，得到最优解。
【算法】分治算法以及用分治算法解决汉诺塔问题
2302_78914800的博客
08-04 248
分治算法的思想是 “ 分而治之 ” ，就是把一个复杂的问题分成两个或更多的相同或相似的子问题，再把子问题分成更小的子问题......知道最后子问题可以简单的直接求解，原问题的解即子问题的解的合并。这个技巧是很多高效算法的基础。
学习笔记--算法（双指针）5
weixin_71085457的博客
08-04 656
学习自用！！
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值