一、时齐Markov链模型

假定 Markov 链是时齐的，也就是转移规律与时间无关，只和时间长短有关，那 Markov 链的状态完全由初始状态和转移概率决定，因此，建立时齐 Markov 链的关键是构造转移概率矩阵，将转移概率矩阵计为:
$$P=\left[\begin{array}{cccc}{p}_{11}& p_{12}& \cdots & p_{1n}\\ p_{21}& p_{22}& \cdots & p_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ p_{n1}& p_{n2}& \cdots & p_{nn}\end{array}\right]$$
已知初始分布${\pi }_0$和转移概率$P$，即可计算各单位时间$n$后的状态分布${\pi }_n$，${\pi }_n={\pi }_0\ast P^n$。

二、案例问题描述

假设2010年全国第六次人口普查中的基线状态数据为$A$，2015年全国人口抽查调查状态数据为$B$，即$B=A\ast P^5$。2010年和2015年人口状态概率$A=(0.4382,0.3933,0.1390,0.0295)$和$B=(0.4051,0.4185,0.1505,0.0260)$。如何求解Markov链模型转移概率矩阵$P$？

三、算法思路&步骤

为了求解$B=A\ast P^5$，可以将矩阵方程转为带约束条件的非线性优化问题。具体如下：
$$\begin{array}{l}\min F=||B-A{P}^5||\\ s.t.\{\begin{array}{l}\sum _{j=1}^{j=4}{p}_{ij}=1,\\ 0<p_{ij}<1,\\ 0<i<4,0<j<4\end{array}\end{array}$$

四、Matlab代码

% 部分Matlab代码---------------------------
clear,clc,close all
%% 计算状态转移矩阵
A = [0.4382,0.3933,0.1390,0.0295];
B = [0.4051,0.4185,0.1505,0.0260];
% 非线性优化fmincon
Aeq = kron(eye(4),ones(1,4)); % 等式约束，行和为1
beq = ones(4,1);
lb = zeros(1,16); % 变量上下界
ub = ones(1,16);
x0 = ones(1,16)*0.1; % 初始解guess
[x, fval] = fmincon(@fun,x0,[],[],Aeq,beq,lb,ub);
x = reshape(x,[],4)';
disp('转移概率矩阵：');disp(x)
%% 计算未来年份
year = 15:5:40;
n = length(year);
data = zeros(n,4);
for i = 1:n
    data(i,:) = A*x^year(i);
end

🍗完整代码传送门