Light bulbs: (Exponential) life time. A light bulb manufacturer wants to estimate the life time of new bulbs. He knows that the life of bulbs follows an exponential distribution. The life times of 10 light bulbs in hours, obtained through an accelerated life test, are 26293, 10123, 3168, 23340, 5459, 13143, 10270, 1699, 15061, 29010. Using a non-informative (improper) flat prior distribution: a) Write the likelihood function. b) Write the posterior distribution. c) Draw the posterior distribution and the likelihood function in the same graph. The manufacturer of these light bulbs wants to specify their guarantee time. This guarantee time should be subject to the next condition: the probability that the time to failure of a light bulb were higher than the guarantee time must be 95%. d) Calculate this guarantee time. \

Now, our knowledge has been enough so far that we have understood that the gamma distribution is a conjugate prior for the exponential distribution. So, let us proceed with our knowledge that the lifetime of the bulb is modeled as an exponential random variable with unknown parameter . \

lifetime_bulbs <- c(26293, 10123, 3168, 23340, 5459, 13143, 10270, 1699, 15061, 29010)
data <- lifetime_bulbs/(24*365)
mean_lifetime <- mean(data)

We will assume that we have no conccrete information about the behaviour of our prior. So, we will choose a non-informative prior. In this case, I suggest the use of

\[ Gamma(\alpha=0.001, \beta = 0.002 ) \] as our choice of prior so that the mean will be $ = = 0.5 $. This will reflect our lack of information we have about the values that can take. If we plot our prior, it looks as follows-

alpha <- 0.01
beta <- 0.02
curve(dgamma(x, shape = alpha, rate = beta), col="orange", lwd=2, lty=2)

Also, if we take a look at the data provided, it becomes a bit tricky to plot and hence I have scaled the data on the scale of years for the lifetime which means, we will be dividing the original values by a factor of (24) x (365). Now, if we consider the data that we have, we can easily towards computing the likelihood using the formuala- \[ p(x| \lambda) = \lambda^{10} exp \{ -\lambda \sum_{i=1}^{10} X_i \} \] and since we assume that our data follows the exponential distribution with parameter \(\lambda\), we can plot the likelihood as follows-

plot( function(x) return(x^10 * exp(-x*sum(data))*2000000), col="darkgreen", lwd=2,
      ylab = "")

Note : We have multiplied our likelihood function with an arbitrary large factor to simplify the plot.

Now, given the likelihood, our posterior function takes the following form- \[ p(\lambda | x) = p(x| \lambda) p(\lambda) \] which can be simplified to, \[ p(\lambda | x) = \lambda^{10} exp\{ -\lambda \sum_{i=1}^{10} X_i \} . Gamma(\alpha=0.001, \beta = 0.002) \]

plot( function(x) return(x^10 * exp(-x*sum(data))*2000000), ylim=c(0,3), 
      lty=2, lwd=2, col="tomato", ylab="", main="Posteror versus Likelihood")
curve(dgamma(x, shape = alpha +length(data), rate = beta + sum(data)  ), 
      col = "steelblue" , lwd=3, add = TRUE, ylab="")
legend("topright", c("posterior","likelihood"),lty=c(1,2), col = c("steelblue", "tomato"),
       lwd=c(2,2))

Now, as we know that our posterior distribution follows the gamma distribution as we have used a non-informative but conjugate prior, we can compute the cdf of failure times using the pgamma function. This will tell us about what value of X is such that- \[ P(X <0.95) = 0.9276907 \]

Similarily, Quantile function qgamma tells that which value of \(\lambda\) is there that makes 95% of the time guarantee that time to failure is greater than that value

qgamma(0.95, alpha +length(data), rate = beta + sum(data))
## [1] 0.9996164