In a nearby forest we can still find some Iberian Lynxs. There are very few of them but we’re not sure how many.

Assuming that \(X\) is the number of lynxs can vary between \(0\) and \(5\) with

\[p(X=x) = {5 \choose x} 0.6^x 0.4^{5-x}\]

for \(x=0,1,\ldots,5\).

Unfortunately, it is very hard to see a wild lynx. Assume that the number of observed lynxs is given by \(Y\) with the following likelihood:

\[p(Y=y|X=x) = \left\{ \begin{array}{ll} {x \choose y} 0.3^y 0.7^{x-y} & 0 \leq y \leq x \\ 0 & x < y \end{array} \right. \]

What would be the conditional distribution of \(X\) after seeing \(2\) lynxs?

linces <- 0:5
priors <- choose(5,linces) * 0.6^linces * 0.4^(5-linces)

likelihood <- function(linces.obs,linces.prior) {
  result = numeric(length(linces.prior))
  for(i in 1:length(linces.prior)) {
    if (linces.prior[i]<linces.obs)
      result[i] <- 0
    else
      result[i] <- (choose(linces.prior[i],linces.obs)
                    *0.3^linces.obs
                    *0.7^(linces.prior[i]-linces.obs))
  }
  return(result)
}

likelihoods <- likelihood(2,linces)

posteriors <- priors * likelihoods / sum(priors * likelihoods)

priors
## [1] 0.01024 0.07680 0.23040 0.34560 0.25920 0.07776
likelihoods
## [1] 0.0000 0.0000 0.0900 0.1890 0.2646 0.3087
posteriors
## [1] 0.0000000 0.0000000 0.1160749 0.3656360 0.3839178 0.1343712
par(mfrow=c(2,1))
barplot(priors,xlab="number of lynxs",
        ylim=c(0,.5),names.arg=linces, 
        main="Prior Distribution")
barplot(posteriors,xlab="number of lynxs",
        ylim=c(0,.5),names.arg=linces, 
        main="After seeing two lynxs")

But giving zero probability in the likelihood can be too strict. What if some error was made when communicating the observations? After the zero, no amount of information will recover the probabilities of having zero or one lynxs.

We could attribute a small probability to consider those rare problems:

\[p(Y=y|X=x) = \left\{ \begin{array}{ll} {x \choose y} 0.3^y 0.7^{x-y} & 0 \leq y \leq x \\ 0.001 & x < y \end{array} \right. \]

linces <- 0:5
priors <- choose(5,linces) * 0.6^linces * 0.4^(5-linces)

likelihood.v2 <- function(linces.obs,linces.prior) {
  result = numeric(length(linces.prior))
  for(i in 1:length(linces.prior)) {
    if (linces.prior[i]<linces.obs)
      result[i] <- 0.001
    else
      result[i] <- (choose(linces.prior[i],linces.obs)
                    *0.3^linces.obs
                    *0.7^(linces.prior[i]-linces.obs))
  }
  return(result)
}

likelihoods <- likelihood.v2(2,linces)

posteriors <- priors * likelihoods / sum(priors * likelihoods)

priors
## [1] 0.01024 0.07680 0.23040 0.34560 0.25920 0.07776
likelihoods
## [1] 0.0010 0.0010 0.0900 0.1890 0.2646 0.3087
posteriors
## [1] 5.729304e-05 4.296978e-04 1.160184e-01 3.654580e-01 3.837309e-01
## [6] 1.343058e-01
par(mfrow=c(2,1))
barplot(priors,xlab="number of lynxs",
        ylim=c(0,.5),names.arg=linces, 
        main="Prior Distribution")
barplot(posteriors,xlab="number of lynxs",
        ylim=c(0,.5),names.arg=linces, 
        main="After seeing two lynxs")

In this case, we see there’s no relevant difference for the remaining \(X=x\) cases.

But let’s assume now that the next 100 observations didn’t see any lynx.

linces <- 0:5
priors <- choose(5,linces) * 0.6^linces * 0.4^(5-linces)

observations <- c(2,rep(0,100))
for(i in 1:length(observations)) {
  likelihoods <- likelihood(observations[i],linces)
  posteriors <- priors * likelihoods / sum(priors * likelihoods)
  priors <- posteriors
}

posteriors
## [1] 0.000000e+00 0.000000e+00 1.000000e+00 1.018860e-15 3.460253e-31
## [6] 3.917237e-47
barplot(posteriors,xlab="number of lynxs",
        ylim=c(0,1),names.arg=linces, 
        main="After 100 no lynx observations")

With the initial likelihood there’s nothing we can do to recover the old hypothesis of having zero or one lynxs. But if we just assign 0.001 instead of 0, the hypothesis of zero lynxs comes to life!

linces <- 0:5
priors <- choose(5,linces) * 0.6^linces * 0.4^(5-linces)

observations <- c(2,rep(0,100))
for(i in 1:length(observations)) {
  likelihoods <- likelihood.v2(observations[i],linces)
  posteriors <- priors * likelihoods / sum(priors * likelihoods)
  priors <- posteriors
}

posteriors
## [1] 1.000000e+00 2.425857e-15 2.118522e-28 2.158478e-43 7.330623e-59
## [6] 8.298755e-75
barplot(posteriors,xlab="number of lynxs",
        ylim=c(0,1),names.arg=linces, 
        main="After 100 no lynx observations")