diff --git a/chapter13/HMM/hmmRecSmoother_.m b/chapter13/HMM/hmmRecSmoother_.m
new file mode 100644
index 0000000..f29148c
--- /dev/null
+++ b/chapter13/HMM/hmmRecSmoother_.m
@@ -0,0 +1,23 @@
+function [ gamma, c ] = hmmRecSmoother_( M, A, s )
+% Forward-backward (recursive gamma no alpha-beta) alogrithm for HMM to compute posterior p(z_i|x)
+% Input:
+%   x: 1xn observation
+%   s: kx1 starting probability of p(z_1|s)
+%   A: kxk transition probability
+%   E: kxd emission probability
+% Output:
+%   gamma: 1xn posterier p(z_i|x)
+%   llh: loglikelihood or evidence lnp(x)
+% Written by Mo Chen sth4nth@gmail.com
+[K,T] = size(M);
+At = A';
+c = zeros(1,T); % normalization constant
+gamma = zeros(K,T);
+[gamma(:,1),c(1)] = nml(s.*M(:,1),1);
+for t = 2:T
+    [gamma(:,t),c(t)] = nml((At*gamma(:,t-1)).*M(:,t),1);  % 13.59
+end
+for t = T-1:-1:1
+    gamma(:,t) = nml(bsxfun(@times,A,gamma(:,t)),1)*gamma(:,t+1);
+end
+