Расширенное уравнение стрижки (цепь Маркова) Python Реализация - PullRequest
Новая
       44

Расширенное уравнение стрижки (цепь Маркова) Python Реализация

0 голосов
/ 03 сентября 2018

Работает над следующей задачей минимизации модели линейного программирования:

enter image description here

Где m - состояние, K - порядок / шаг, X - матрица вероятностей устойчивого состояния.

Мне нужно выяснить значения лямбда-вектора при минимизации суммы значений в векторе W .

Ниже приведен пример и код, который я пытаюсь использовать: 

import numpy as np
from pulp import *

state_sample = [1,1,2,2,1,3,2,1,2,3,1,2,3,1,2,3,1,2,1,2]

def transition_matrix(transitions,order):
    n = len(state) #number of states

    M = [[0]*n for _ in range(n)]

    for (i,j) in zip(transitions,transitions[order:]):
        M[i-1][j-1] += 1

    #now convert to probabilities:
    for row in M:
        s = sum(row)
        if s > 0:
            row[:] = [f/s for f in row]
    return M

One_Step_Transition_Matrix= transition_matrix(state_sample,order=1)
one_step_array = np.array(One_Step_Transition_Matrix)

two_step_array = np.array(transition_matrix(state_sample,order = 2))

def steady_state_prop(p):
    dim = p.shape[0]
    q = (p-np.eye(dim))
    ones = np.ones(dim)
    q = np.c_[q,ones]
    QTQ = np.dot(q, q.T)
    bQT = np.ones(dim)
    return np.linalg.solve(QTQ,bQT)

steady_state = steady_state_prop(one_step_array)

Q_Arr = np.vstack((np.matmul(one_step_array,steady_state),np.matmul(two_step_array,steady_state))).transpose()

Weight_vec = []
Number_of_states = Q_Arr.shape[0]
for x in range(Number_of_states):
    Weight_vec.append('w'+str(x+1))

L1 = LpVariable("L1",0,None)
L2 = LpVariable("L2",0,None)

prob = LpProblem("Problem",LpMinimize)

for s in range(Number_of_states):
    Weight_vec[s] = LpVariable('w'+str(s+1),0,None)
count = 0

for row in Q_Arr:
    prob += Weight_vec[count] >= steady_state[count] - row[0]*L1 - row[1]*L2
    prob += Weight_vec[count] >= -steady_state[count] + row[0]*L1 + row[1]*L2
    count = count + 1

prob += L1 >= 0
prob += L2 >= 0

prob += L1 + L2 == 1

for s in range(Number_of_states):
    prob += Weight_vec[s] >= 0

#objective
prob += sum(Weight_vec)

status = prob.solve(GLPK(msg=0))
LpStatus[status]

result = []

for s in range(Number_of_states):
    result.append(value(Weight_vec[s]))
result.append(value(L1))
result.append(value(L2))

print (result)

Вот результат из приведенного выше кода: 

[0.00854214, 0.083957, 0.167995, 1.0, 0.0]

Если я положу значения W в уравнении выше, Мои уравнения не будут удовлетворены.

Не могли бы вы помочь мне!

...