用条件python填充矩阵

Question

我正在将此 matlab 代码转换为 Python，但显然我犯了一个小错误，即结果不一样。这是matlab代码：

alpha = 0.1;
beta = 0.01;
Dt = 0.001;
N = 1000;
time = 0:Dt:40; 
state = zeros(length(time),N); 
R = rand(length(time)-1,N);

for j = 1:N
    for i = 2:length(time)
        if (state(i-1,j) == 0) 
            if (R(i-1,j) <= alpha*Dt) 
                state(i,j) = 1;
            else
                state(i,j) = 0;
            end;
        else                     
            if (R(i-1,j) <= beta*Dt) 
                state(i,j) = 0;
            else
                state(i,j) = 1;
            end; 
        end;
    end; 
end;
Open_stochastic = sum(state,2)/N; 
open_deterministic = (1-exp(-time*(alpha+beta)))*alpha/(alpha+beta);

plot(time, Open_stochastic, ’k-’, time, Open_deterministic, ’k--’)

这是我的python代码：

import numpy as np 
import matplotlib.pyplot as plt 
%matplotlib inline 

alpha = 0.1
beta = 0.01
Dt = 0.001
N = 1000

time = np.arange(0, 40, Dt)
state = np.zeros((len(time), N))
R = np.random.rand(len(time), N)

for j in range(0, N-1):
    for i in range(1, len(time)-1):
        if state[i-1][j] == 0:
            if R[i-1][j] <= alpha*Dt:
                state[i][j] = 1
            else:
                state[i][j] = 0
            if R[i-1][j] <= beta*Dt:
                state[i][j] = 0
            else:
                state[i][j] = 1

    
open_stochastic = state.sum(1)/N 
plt.plot(time, open_stochastic)

open_stochastic 元素全为零，绘图只是一条水平线，但在 matlab 代码中，绘图是指数的（参见图片）

Answer 1

下一个问题是您在最里面的if/else 对之间缺少一个else:。此外，您应该从两个for 语句中删除-1。请记住，Python 范围不包括 endpint。 range(N) 已经从 0 运行到 N-1。通过这些变化，我得到了一些接近你想要的东西。

更新

这是正确的代码：

...
for j in range(N):
    print(j,end='\r')
    for i in range(1, len(time)):
        if state[i-1][j] == 0:
            if R[i-1][j] <= alpha*Dt:
                state[i][j] = 1
            else:
                state[i][j] = 0
        else:    # <<< this is the line you were missing.
            if R[i-1][j] <= beta*Dt:
                state[i][j] = 0
            else:
                state[i][j] = 1
    
open_stochastic = state.sum(1)/N 
plt.plot(time, open_stochastic)
plt.show()