Contents

Homework 13

Homework 13

  • Type your name and email in the “Student details” section below.

  • Develop the code and generate the figures you need to solve the problems using this notebook.

  • For the answers that require a mathematical proof or derivation you can either:

    • Type the answer using the built-in latex capabilities. In this case, simply export the notebook as a pdf and upload it on gradescope; or

    • You can print the notebook (after you are done with all the code), write your answers by hand, scan, turn your response to a single pdf, and upload on gradescope.

  • The total homework points are 100. Please note that the problems are not weighed equally.

Note

  • This is due before the beginning of the next lecture.

  • Please match all the pages corresponding to each of the questions when you submit on gradescope.

Student details

  • First Name:

  • Last Name:

  • Email:

import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set(rc={"figure.dpi":100, 'savefig.dpi':300})
sns.set_context('notebook')
sns.set_style("ticks")
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('retina', 'svg')
import numpy as np
import scipy.stats as st

Problem 1 - Maximum likelihood estimate of an Exponential random variable

You are managing toll station and your job is to come up with a probabilistic model for the time that passes between car arrivals. You have at your disposal data:

\[ T_{1:N} = (t_1, t_2, \dots, t_N). \]

Here \(t_1\) is the time that passed until the first event, \(t_2\) is the time that passed from the first event to the second, and so on. You decide to model these random time intervals using an Exponential distribution with uknown rate parameter \(\lambda\), i.e.,

\[ T_i | \lambda \sim \text{Exponential}(\lambda). \]

In terms of the PDF:

\[ p(t_i | \lambda) = \lambda e^{-\lambda t_i}. \]

You have decided to use maximum likelihood to fit the parameter \(\lambda\). Answer the following questions:

  • Find the mathematical form of the likelihood of the data \(p(t_{1:N}|\lambda)\).

Answer:





  • Find the mathematical form of the log-likelihood \(J(\lambda) = \log p(t_{1:N}|\lambda)\).

Answer:





  • Take the derivative \(\frac{dJ(\lambda)}{d\lambda}\) of \(J(\lambda)\) with respect to \(\lambda\).

Answer:





  • Find the MLE \(\hat{\lambda}\) by solving the equation \(\frac{dJ(\lambda)}{d\lambda} = 0\) for \(\lambda\).

Answer:





Problem 2 - Failure of a mechanical component

Assume that you designing a gear for a mechanical system. Under normal operating conditions the gear is expected to fail at a random time. Let \(T\) be a random variable capturing the time the gear fails. What should the probability density of \(T\) look like? Well, when the gear is brand new, the probability density should be close to zero because a new gear does not fail under normal opearating conditions. As time goes by, the probability density should increase because various things start happening to the material, e.g., crack formation, fatigue, etc. Finally, the probability density must again start going to zero as time further increases because nothing lasts forever… A probability distribution that is commonly used to model this situation is the Weibull. We are going to fit some fail time data to a Weibull distribution and then you will have to answer a few questions about failing times.

# Time to fail in years under normal operating conditions
# Each row is a different gear
time_to_fail_data = np.array([
    10.5,
    7.5,
    8.1,
    8.4,
    11.2,
    9.3,
    8.9,
    12.4
])

Here is a Weibull distribution fitted to the data using maximum likelihood:

fitted_params = st.exponweib.fit(time_to_fail_data, loc=0)
T = st.exponweib(*fitted_params)

Let’s plot the probability density of this:

fig, ax = plt.subplots(dpi=100)
ts = np.linspace(0.0, 20.0, 100)
ax.plot(ts, T.pdf(ts))
ax.set_xlabel('$t$ (years)')
ax.set_ylabel('$p(t)$');

  • Find the mean fail time and its variance. Hint: Do not integrate anything by hand. Just use the functionality of scipy.stats.

# Your code here
t_mean = 0.0 # Change me
t_var = 0.0 # Change me
print('E[T] = {0:1.2f}'.format(t_mean))
print('V[T] = {0:1.2f}'.format(t_var))

  • Plot the cumulative distribution function of \(T\).

# Your code here

  • Plot the probability that gear survives for more than \(t\) as a function of \(t\). That is, plot the function:

\[ S(t) = p(T > t). \]

Hint: First express this function in terms of the cumulative distribution function of \(T\).

# Your code here

  • Find the probability that the gear lasts anywhere between 8 and 10 years.

# Your code here

  • If you were to sell the gear, how many years “warranty” would you offer and why?

Answer:





previous

Predictive checking

next

Lecture 14: Covariance, correlation, and linear regression with one variable