Как построить эту интегральную функцию с границей бесконечности?

Question

Я пытаюсь построить функцию Маркума, которая является интегралом от функции Гаусса на границе.

Это функция:

<a href="https://www.codecogs.com/eqnedit.php?latex=Q(y)=\int_{y}^{\infty}\frac{1}{2\pi}\exp(-\frac{x^2}{2})~dx" target="_blank"><img src="https://latex.codecogs.com/gif.latex?Q(y)=\int_{y}^{\infty}\frac{1}{2\pi}\exp(-\frac{x^2}{2})~dx" title="Q(y)=\int_{y}^{\infty}\frac{1}{2\pi}\exp(-\frac{x^2}{2})~dx" /></a>

y - это значение с плавающей точкой, которое пользователь должен ввести. Я написал код для вычисления этой функции при желаемом значении. Однако то, что я собираюсь построить, выглядит примерно так:

! https://ibb.co/SRMsfY3

Вот код для вычисления значения Q (y):

from scipy.integrate import quad
import scipy, numpy, math
from scipy import integrate


def integrand(x):
   Q_y = (1/(math.sqrt(2*math.pi)))*math.exp((-x**2)/2)
   return Q_y


y = input('y=')
ans, err = quad(integrand, float(y), math.inf)
print(ans)

Я пытался построить его, используя следующий код:

x_data=np.arange(-20,20,0.1)    
z_data=np.arange(float(y),1000,1)

    ans, err = quad(integrand, float(y), math.inf)
    print(ans)

    for item in x_data:
        z_data[i]=integrand(x_data[i])[0]
        i+=1

    fig = plt.figure(1, figsize=(4,4))
    plt.plot(x_data,z_data,color= 'blue',linestyle='--', label=label1)

, который возвращает следующую ошибку:

 File "C:/Users/m.rafiee/.spyder-py3/temp.py", line 59, in <module>
    z_data[i]=integrand(x_data[i])[0]

IndexError: index 400 is out of bounds for axis 0 with size 400

Я также не знаю, как добавить аннотации, показанные на рисунке, и область хеширования при построении.

Я очень благодарен за любую помощь заранее.

Answer 1

Векторизация вашей функции и сюжета:

from scipy.integrate import quad
import scipy as sc, numpy as np, math
from scipy import integrate
import matplotlib.pyplot as plt

def integrand(x):
   Q_y = (1/(math.sqrt(2*math.pi)))*math.exp((-x**2)/2)
   return Q_y

func = np.vectorize(integrand)

x = np.arange(-3, 3, 0.1)

fig = plt.figure()
ax = plt.gca()

y = input('y=')

y = float(y)

result = integrand(y)

x_user = np.arange(y, 3, 0.1)

ax.axvline(y, color='green', lw=2, alpha=0.5)

ax.fill_between(x_user, func(x_user), facecolor="none", hatch="/", edgecolor="red")

ax.annotate('Q({0})={1}'.format(y, round(result, 2)), xy=(0, result), xytext=(0, result))

ax.plot(x, func(x))