Как совместить интеграционное уравнение в Python?

Question

Я рассчитываю ключевую скорость (R ^ Rate-wise), интегрируя R (eta) по всем возможным eta от 0 до 1, с распределением вероятностей (PDTC), которое является логарифмически нормальным распределением.

Уравнение логарифмического нормального распределения:

Уравнение R (eta):

Следовательно, R ^ по ставкам = Integrate_0 ^ 1 (R (эта) * P (эта) * д эта):

Это Python-код нормального логарифмического распределения:

x=np.linspace(0,1,1000)
sigma0=[0.9]
color=['green']
for i in range(len(sigma0)):
    sigma=sigma0[i]
    y=1/(x*sigma*np.sqrt(2*np.pi))*np.exp(-(np.log(x/0.3)+(1/2*sigma*sigma))**2/(2*sigma*sigma))
    plt.plot(x,y,color[i])
plt.title('Lognormal distribution')
plt.xlabel('x')
plt.ylabel('lognormal density distribution')
#plt.xlim((0,0.002))
plt.ylim((0,5))
plt.show()

Это Python-код R (eta):

n1=np.arange(10, 55, 1)
n=10**(-n1/10)


Y0=1*(10**-5)
nd=0.25
ed=0.03
nsys=nd*n
QBER=((1/2*Y0)+(ed*nsys))/(Y0+nsys)
H2=-QBER*np.log2(QBER)-(1-QBER)*np.log2(1-QBER)
Rsp=np.log10((Y0+nsys)*(1-(2*H2)))
print (Rsp)

plt.plot(n1,Rsp)
plt.xlabel('Loss (dB)')
plt.ylabel('log10(Rate)')
plt.show()

Мой вопрос: какинтегрировать R (ETA) по возможным ETA от 0 до 1?Выходные данные должны быть на следующем рисунке (R ^ по шкале):

Ссылку на статью можно найти по этой ссылке: https://arxiv.org/pdf/1712.08949.pdf

Answer 1

Я не понимаю, в какой диапазон значений eta вы пытаетесь интегрировать, поэтому я выбрал x = np.linspace(0,2,1000) и немного организовал код, чтобы он подходил для scipy.integrate.quad.

x=np.linspace(0,2,1000)

def log_normal(x,sigma):
    y=1/(x*sigma*np.sqrt(2*np.pi))*np.exp(-(np.log(x/0.3)+(1/2*sigma*sigma))**2/(2*sigma*sigma))

    return y


def R(x,nd,Y0,ed):
    nsys = x*nd
    QBER=((1/2*Y0)+(ed*nsys))/(Y0+nsys)
    H2=-QBER*np.log2(QBER)-(1-QBER)*np.log2(1-QBER)
    out = (Y0+nsys)*(1-(2*H2))
    return out


def out(x,nd,Y0,ed,sigma):

    term1 = log_normal(x,sigma)
    term2 = R(x,nd,Y0,ed)

    return term1*term2

hyp = (0.25,1*(10**-5), 0.03,0.9) #nd,Y0,ed,sigma


lg = log_normal(x[1:],hyp[3])

r = R(x[1:],hyp[0],hyp[1],hyp[2])
plt.plot(lg,label = "log_normal sigma = 0.9")
plt.plot(r,label = "R")
plt.legend()
plt.show()

plt.plot(np.log10 ( out(x,hyp[0],hyp[1],hyp[2],0.9) ),label = "sigma = 0.9" )
plt.plot(np.log10 ( out(x,hyp[0],hyp[1],hyp[2],0.3) ) ,label = "sigma = 0.3")
plt.plot(np.log10 ( out(x,hyp[0],hyp[1],hyp[2],0.6) ) ,label = "sigma = 0.6")

#plt.xlim((0,50))
#plt.ylim((-10,0))
plt.title("Rate")
plt.legend()
plt.show()

И теперь мы можем интегрировать.

from scipy.integrate import quad

sol = quad(out,a = 0,b =1 ,args=hyp)
sol
#(0.03721392434457473, 1.9175746760101603e-10)