In [1]:
import numpy as np
import scipy.optimize
import matplotlib.pyplot as plt
In [2]:
f = lambda a: np.exp(-24*a) - 0.5 * np.exp(-50*a) - 0.5
In [3]:
fig, ax = plt.subplots(figsize=(8,5))
a_min, a_max = 0, 0.01
a = np.linspace(0, 0.01, 100)
ax.plot(a, f(a))
ax.hlines(0, a_min, a_max, linestyles=':')
plt.show()
In [4]:
a_star = scipy.optimize.brentq(f, 0.002, 0.004)
a_star
Out[4]:
0.0032034184388007483
In [5]:
CE_a = 20
In [6]:
def CE(x, p, a):
    x = np.asarray(x)
    p = np.asarray(p)
    return -np.log(p @ np.exp(-a * x)) / a
In [7]:
CE([20], [1], a_star)
Out[7]:
19.999999999999982
In [8]:
CE_b = CE([50, -10], [0.6, 0.4], a_star)
CE_b
Out[8]:
24.600322886059406
In [9]:
CE_c = CE([50, 0, 10], [0.4, 0.3, 0.3], a_star)
CE_c
Out[9]:
22.204239619220147
In [ ]: