模拟非马尔可夫随机过程

2024-05-14 14:04:29 发布

您现在位置:Python中文网/ 问答频道 /正文

考虑一个非马尔可夫过程,它描述了从状态a到状态B的转变,由一个连续时间内的响应时间分布控制,并有额外的恒定流入状态a和从状态B流出。参见这里的过程说明。See process illustration here.

例如,我有100个A型单元格,它们转换为B型,其中事件之间的分布是非指数分布的。为了模拟这种情况,我需要逐个查看每个单元格,并检查何时发生状态转换。我努力将恒定的通量整合到状态A中,同时也考虑到状态B的转变,我是否需要单独关注每个细胞,或者我是否可以从指数分布中得出另一个细胞何时加入系统?我的脑子完全卡住了。你知道吗

def semi_markov(start, stop, nsteps, ncells, nstates = 3):

# initialize cells
cells = np.zeros((ncells, nsteps, nstates))    
time = np.linspace(start, stop, nsteps)

# generate random number for each cell to compare with integral
# of response-time distribution
numbers_rnd = [np.random.rand() for i in range(ncells)]
fate_times = [np.random.exponential(1.) for i in range(ncells)]

# for each time point loop over each cell and check if transition occurs
for i in range(len(time)-1):
    t = time[i]
    t_new = time[i+1]

    cell_j = cells[:,i,:]

    for j in range(cells.shape[0]):          
        cell = cell_j[j,:]           

        n_rnd = numbers_rnd[j]
        fate_time = fate_times[j]

        # if cell is type A check if transition to B occurs
        if cell[0] == 0:
            if n_rnd > cell[1]:
                cell[1] = cell[1]+ (non_exp_cdf(t_new)-non_exp_cdf(t))
             else:
                cell[0] = 1
                cell[2] = t

        # if cell is type B check if transition occurs    
        if cell[0] == 1 and fate_time < (t - cell[2]):
            cell[0] = 2                

        cells[j,i+1,:] = cell

return [cells, time]

Tags: inforiftime状态npcellrange

热门问题