如何模拟随机过程，直到路径的所有元素都是正数?

Question

我想使用 Python 在连续的时间内模拟某个随机过程的路径。我写了以下函数来模拟它(当然np指的是numpy)：

def simulate_V(v0, psi, theta, dt, xsi, sample, diff = False):
_, dB = independent_brownian(dt, 1, sample, diff = True)
path_lenght = len(dB) + 1
V = np.zeros(path_lenght)
dV = np.zeros(len(dB))
V[0] = v0
for i in range(len(dV)):
dV[i] = psi * (theta - V[i]) * dt + xsi * np.sqrt(V[i]) * dB[i]
V[i+1] = V[i] + dV[i]

if diff == True:
return(V, dV)
else:
return(V)

independent_brownian函数只是为标准布朗运动创建一条路径。为了完整起见，这里是：

def independent_brownian(dt, n, sample, diff=False):
'''
Creates paths of independent Brownian motions. Returns increments if diff == True.
--------
dt -> variance of increments\
n -> how many paths\
sample -> length of sampled path\
diff -> return increments or not?\
Returns a (n, sample)-shaped array with the paths
'''
increments = np.sqrt(dt) * np.random.randn(n, sample)
paths = np.cumsum(increments, axis=1)
b = np.zeros((n, sample + 1))
b[:, 1:] = paths
if n==1:
b = b.flatten()
increments = increments.flatten()
if diff == True:
return(b, increments)
else:
return(b)

碰巧我的模型背后的数学意味着过程 $V_t$(由上面的代码中的离散化V表示)必须是正数。但是，数组dB可能包含绝对值较大的负元素的情况。我想自动执行以下"选择"程序：

1. 尝试遵循simulate_V中描述的循环;
2. 在某些时候，如果V[i]低于零，中断该过程，dB对另一个序列进行采样并重新开始;
3. 当V的所有元素都是正的时停止;

自动化此过程的好方法是什么？现在，我明白，如果V[i]低于零，我将在 numpy 中得到一个nan，但它不会抛出错误或停止该过程。

Answer 1

为了概括起见，并且由于我对布朗运动不是很熟悉，我将为我们有变量dB的情况实现一个通用机制，我们用它来创建另一个变量V，我们需要重复这一代V，直到某个条件成立。

import numpy as np
dB = np.random.normal(size=3) # initializing dB
found_sufficient_v_flag = False # we define a flag that would tell us when to stop
while not found_sufficient_v_flag: # while this flag is False, we continue searching
v = np.zeros(3) # initializing v
for i in range(3): # for loop
v[i] = dB[i] + np.random.normal() # we change the v[i] element
if v[i] < 0: # if it is negative, there is no point in continuing, so we break the loop
dB = np.random.normal(size=3) # but we need to instantiate a different dB
break 
if np.all(v > 0): # we arrive here if either v[i] < 0, or all v[i] > 0. in the latter case we stop
found_sufficient_v_flag = True
print(dB)
print(v)

这为我提供了以下输出：

[2.27582634 0.77008881 0.28388536]
[2.55101104 3.10944337 0.55829105]

你可以看到V的条件确实成立。

如果需要任何其他澄清，请告诉我。