Markov Chain

http://www.52nlp.cn/lda-math-mcmc-%E5%92%8C-gibbs-sampling1

Java：

void markovConverge() {
        double x[][] = new double[][] { { 0.6500, 0.2800, 0.0700 }, { 0.1500, 0.6700, 0.1800 },
                { 0.1200, 0.3600, 0.5200 } };

        // todo:使用matlab来写
        final int ITERARTION = 100000;
        // double pi[] = new double[] { 0.7, 0.2, 0.1 };
        // double sample[] = new double[10000];
        int walkStates[] = new int[ITERARTION];
        int pre_state = 2;// top
        // sample[0] = pi[pre_state];
        for (int i = 1; i < ITERARTION; i++) {
            double u = Math.random();// 换成norm呢？
            int state;
            if (u < x[pre_state][0]) {
                state = 0;// top
            } else if (u < x[pre_state][1] + x[pre_state][0]) {
                state = 1;// mid
            } else {
                state = 2;// low
            }
            // sample[i] = pi[pre_state] * x[pre_state][state];
            walkStates[i] = pre_state;
            pre_state = state;
        }

        double p1, p2, p3;
        p1 = p2 = p3 = 0;
        for (int i = 0; i < ITERARTION - 1; i++) {
            if (walkStates[i] == 0) {
                p1 += 1.0 / (ITERARTION - 1);
            } else if (walkStates[i] == 1) {
                p2 += 1.0 / (ITERARTION - 1);
            } else {
                p3 += 1.0 / (ITERARTION - 1);
            }
        }
        System.err.println(p1 + "  " + p2 + "  " + p3);
    }

我写了段matlab代码，通过Q，求出pi。

Q = [0.65 0.28 0.07;0.15 0.67 0.18;0.12 0.36 0.52];
ITERARTION = 100000;
walkStates=zeros(ITERARTION,1);
pre_state = 2;

for i = 1:ITERARTION
u = rand;
if u <= Q(pre_state,1)
state = 1;
elseif u <= Q(pre_state,1) + Q(pre_state,2)
state = 2;
else
state = 3;
end
walkStates(i) = pre_state;
pre_state = state;
end

begin = 90000;

pi = [length(find(walkStates(begin:ITERARTION,:)==1))/(ITERARTION-begin),length(find(walkStates(begin:ITERARTION,:)==2))/(ITERARTION-begin),length(find(walkStates(begin:ITERARTION,:)==3))/(ITERARTION-begin)]