install.packages('Ryacas')
  1. A bin of 1000 turnbuckles has an unknown number D of defectives. A sample of 100 turnbuckles has 2 defectives. The maximum likelihood estimate for D is the number of defectives which gives the highest probability for obtaining the number of defectives observed in the sample. Guess this number D and then write a computer program to verify your guess.

Solution:

For a single turnbuckles, Let X be a random variable with value 1 if the turnbuckles is defective and 0 if turnbuckles is nondefective. The probability function of X:

\(h(N, k, n, x)=\frac{\binom{k}{x}\binom{N-k}{n-x}}{\binom{N}{n}} (0 \leq x \leq n)\)

As $= $

So \(h_k > h_{k-1}\) if \(k<x(N+1)/n\) So \(h_k < h_{k-1}\) if \(k>x(N+1)/n\)

Thus, \(h_k\) takes its maximum value when its maximum value when k is the largest interger\(\frac{x}{n}(N+1)\). In 1000 turnbuckles, a sample of 100 turnbuckles has 2 defectives. The maximum likelihood estimate of the total number of defective, number D, in the sample is 20 or 2%.

In a size n sample, the probability of finding D defectives,\({\Pi}\), is:
\(p = \frac{n!}{D!(n-D)!}\Pi^D(1-\Pi)^{n-D}\)

\(\frac{\delta p}{\delta \hat{{\Pi}}}=B(\frac{\hat{n}!}{D!(n-D)!}\hat{\Pi}^{D-1}(1-\hat{\Pi})^{n-D})-(n-D)(\frac{\hat{n}!}{D!(n-D)!}\hat{\Pi}^D(1-\hat{\Pi})^{n-D-1})=0\)

\(D\hat{\Pi}^{D-1}(1-\hat{\Pi})^{n-D}-(n-D)\hat{\Pi}^D(1-\hat{\Pi})^{n-D-1}=0\)

\(D\hat{\Pi}^{-1}=(n-D)\hat{\Pi}^D(1-\hat{\Pi})^{-1}\)

\(D\hat{\Pi}^{-1}=(n-D)(1-\hat{\Pi})^{-1}\)

\(D(1-\hat{\Pi})^{-1}=(n-D)\hat{\Pi}^{-1}\)

\(\hat{\Pi}=D/n\)

library(Ryacas) 
p <- Sym("p") 
s <- expression(factorial(n)/(factorial(D)*factorial(n-D))*p^D*(1-p)^(n-D)) 
deriv(s,p)
## expression({
##     .expr3 <- n - D
##     .expr6 <- factorial(n)/(factorial(D) * factorial(.expr3))
##     .expr8 <- .expr6 * p^D
##     .expr9 <- 1 - p
##     .expr10 <- .expr9^.expr3
##     .value <- .expr8 * .expr10
##     .grad <- array(0, c(length(.value), 1L), list(NULL, c("p")))
##     .grad[, "p"] <- .expr6 * (p^(D - 1) * D) * .expr10 - .expr8 * 
##         (.expr9^(.expr3 - 1) * .expr3)
##     attr(.value, "gradient") <- .grad
##     .value
## })

So we know

\(\frac{n!}{D!(n-D)!}D*p^{D-1}(1-p)^{n-D}-\frac{n!}{D!(n-D)!}(n-D)p^D((1-p)^{n-D}-1)=0\)

\(D*p^{D-1}(1-p)^{n-D}-(n-D)p^D((1-p)^{n-D}-1)=0\)

\(p=D/n\)