<jupyter_start><jupyter_text>Problem 3  Question a)Let us write the likelihood function in this setting : $$L(\alpha, \lambda) = f_{\alpha, \lambda}(T_{1}) \times f_{\alpha, \lambda}(T_{2}) \times \mathbb{P}(T > 100K)\times \mathbb{P}(T > T_4)\times \mathbb{P}(T > T_5)$$After calculations, the negative log-likelihood is equal to : $$-l(\alpha, \lambda) = -[2\log(\alpha \lambda) + (\alpha - 1)\log(\lambda^{2}T_{1}T_{2}) - (\lambda)^{\alpha}(T_{1}^{\alpha} + T_{2}^{\alpha})] + (\lambda)^{\alpha}(T_{3}^{\alpha}+T_{4}^{\alpha} + T_{5}^{\alpha})$$where $T_{1} = 44K, T_{2} = 26K, T_{3} = 100 K, T_{4} = 19K, T_{5} = 45K$. We can solve the MLE numerically with `scipy` module `optimize`.<jupyter_code>import numpy as np
from scipy import optimize
from scipy.stats import weibull_min
from math import gamma

# Problem 3
T1 = 44 * 10 ** 3
T2 = 26 * 10 ** 3
T3 = 100 * 10 ** 3
T4 = 19 * 10 ** 3
T5 = 45 * 10 ** 3


def neg_log_L(gamma, alpha) -> float:
    A = 2 * np.log(alpha / gamma) + (alpha - 1) * np.log(T1 * T2 / (gamma ** 2)) - \
        ((1 / gamma) ** alpha) * ((T1 ** alpha) + (T2 ** alpha))
    B = -((1 / gamma) ** alpha) * (T3 ** alpha + T4 ** alpha + T5 ** alpha)
    return -A - B


## Numerically solve MLE equations.
epsilon = 1 * 10 ** -12
x0 = [10, 10]
bnds = ((epsilon, np.inf), (epsilon, np.inf))
fun = lambda x: neg_log_L(x[0], x[1])
solver = optimize.minimize(fun, x0=x0, bounds=bnds)
print(" {} iterations \n".format(solver.nit),
      "lambda = {} and alpha = {}".format(1 / solver.x[0], solver.x[1]))<jupyter_output>171 iterations 
 lambda = 1.0823481787828909e-05 and alpha = 1.535248714292919<jupyter_text>Question b) We want to compute a $(1-\delta)\%$ confidence interval for  $\mathbb{E}(T)$. Theoretically : $$ \mathbb{E}(T) = \frac{1}{\lambda} \Gamma(1+ \frac{1}{\alpha})$$So we have the plug-in estimate of $\mathbb{E}(T)$ by using $\hat \lambda_{MLE}$ and $\hat \alpha_{MLE}$. Let us compute a parametric bootstrap confidence interval. - 1) Generate $(\mu_{n}^{(1)}, \ldots, \mu_{n}^{(m)} )$ where $\mu_{n}^{(i)} = \frac{1}{n} \displaystyle \sum_{k=1}^{n}T_{k}^{(i)}$- 2) Find $x$ and $y$ such that $\mathbb{P}(\mu_{n} - \mu_{*} < x ) = 1 - \delta/2$ and $\mathbb{P}(\mu_{n} - \mu_{*} < y ) = \delta/2$       where $\mu_{*}$ is the plug-in estimate obtained with the MLE.<jupyter_code>## Parametric Boostrap (1-delta) confidence interval
lmda = 1 / solver.x[0]
alpha = solver.x[1]
# reference parameter
mu_star = (1 / lmda) * gamma(1 + 1 / alpha)
print(" In average, {} kms before a reliability problem occurs (plug-in estimation with MLE)".format(mu_star))

m = 100  # number of iterations to aggregate
bootstrap_estimates = []
for i in range(m):
    T_bootstrap = weibull_min.rvs(c=alpha, scale=1 / lmda, size=100)
    bootstrap_estimates.append(np.mean(T_bootstrap))

delta = 0.1
print(" delta = {}".format(delta), "\n", "m = {} ; ".format(m), "n = 100")
# Upper bound : we find x s.t prob(estimate - mu_star < x) = 1-delta/2
x = 2000
count = 0
while (count / m != 1 - delta / 2):
    count = 0
    for i in range(m):
        if bootstrap_estimates[i] - mu_star < x:
            count += 1
    # print(count / m)
    x += 10
print(" Prob(mu_n - mu_estimate < {}) = {}".format(x, count / m))

# Lower bound : we find y s.t prob(estimate - mu_star < y ) = delta/2
y = -10000
count = 0
while (count / m != delta / 2):
    count = 0
    for i in range(m):
        if bootstrap_estimates[i] - mu_star < y:
            count += 1
    # print(count / m)
    y += 10
print(" Prob(mu_n - mu_estimate < {}) = {}".format(y, count / m))

print(" Parametric Boostrap {}% confidence interval : [{},{}]".format((1 - delta)*100, -x + np.mean(bootstrap_estimates),
                                                                      - y + np.mean(bootstrap_estimates)))<jupyter_output>In average, 83182.41571001986 kms before a reliability problem occurs (plug-in estimation with MLE)
 delta = 0.1 
 m = 100 ;  n = 100
 Prob(mu_n - mu_estimate < 8170) = 0.95
 Prob(mu_n - mu_estimate < -9780) = 0.05
 Parametric Boostrap 90.0% confidence interval : [74749.82839853111,92699.82839853111]<jupyter_text>Problem 6Written part is in the written section of the report.<jupyter_code># Problem 4
## Posterior distribution
X = np.array([6.00, 4.82, 3.35, 2.38, 3.59, 4.12, 4.98, 2.69, 6.24, 6.77,
              6.22, 5.42, 5.42, 3.10, 4.65, 4.24, 4.53, 4.62, 5.36, 2.57])

Sigma = (1 + (0.5) * np.sum(X[0:19]**2))**-1
X_i_1 = X[0:19]
X_i = X[1:20]
Mu = Sigma * ( (0.5) * np.sum(X_i_1*(X_i-2)) + 0.5)
print("rho*|X follows a normal distribution with mean = {} and variance = {}".format(Mu, Sigma))

## Prob(X_23 > 4)
n = 100000
X23_estimates = []
count = 0
for i in range(n) :
    rho = np.random.normal(Mu, Sigma, 1)
    epsilon = np.random.normal(0,2,3)
    X23 = (rho**3) * X[19] + 2 * (1 + rho + rho**2) + epsilon[0] + epsilon[1] * rho + epsilon[2] * rho**2
    X23_estimates.append(X23)
    if X23 >= 4 :
        count += 1

print("X23 estimate = {}".format(np.mean(X23_estimates)))
print("Prob(X_23 > 4) = {}".format(count/n))<jupyter_output>rho*|X follows a normal distribution with mean = 0.5160781386089336 and variance = 0.0045245286006764165
X23 estimate = 3.9179892166389694
Prob(X_23 > 4) = 0.4861