Я хотел реализовать формулу, найденную в статье.Вот формула:

(В этой формуле F - это распределение вероятностей, я взял гауссовское для своих тестов).
У меня много проблем:
-Самая важная из них, это то, что моя функция python возвращает отрицательные значения, что не должно происходить, так как формула является вероятностью, поэтому она должна быть положительной
-Я получаю ошибку от numpy, говоря:
Обнаружена ошибка округления, которая препятствует достижению требуемого допуска.Ошибка может быть недооценена.return 1- (mu / lambdaf) ** (k / 2) * quad (tointegrate, 0.0001, t) [0] C: /Users/NLucis/.spyder-py3/temp.py: 32: IntegrationWarning: максимальное числоподразделений (50) был достигнут.Если увеличение лимита не приводит к улучшению, рекомендуется проанализировать подынтегральное выражение, чтобы определить трудности.Если можно определить положение локальной сложности (особенность, разрыв), то, вероятно, выиграет от разбиения интервала и вызова интегратора на поддиапазонах.Возможно, следует использовать специальный интегратор.
- Захват странный (слишком крутой, возможно, указывает на ошибку в реализации)

Так что я, должно быть, сделал ошибку в реализации, какая?
def survival_time(k,f,t, lambda0, mu, F) :
lambdaf = lambda0*F(f)
tointegrate = (lambda u : math.exp((lambdaf-mu)*u) * (k/u) * scipy.special.iv(k, 2*math.sqrt(mu*lambdaf)*u) )
return 1-(mu/lambdaf)**(k/2) * quad(tointegrate,0.0001,t )[0]
def F(z) :
return 0.5*(1 + scipy.special.erf(z/math.sqrt(2)))
# setting the x - coordinates
x = np.arange(0, 1000, 10)
# setting the corresponding y - coordinates
vst = np.vectorize(survival_time)
y = vst(10, 10, x, 1, 1, F)
# potting the points
plt.plot(x, y)
# function to show the plot
plt.show()