Reinforcing Likelihood concepts in L2.3-2.6 and developing R skills.

Problem 1: Sampling Distribution vs. Likelihood function In an experiment, 10 seed are planted. 7 are observed to germinate and produce seedlings.

i. What is an appropriate distribution to use to model this process? Why?

A binomial would be an appropriate distribution because we expect only two outcomes (germinate, fail), they are mutually exclusive and the results are independent.

ii. Plot the Sampling distribution: Calculate P(x|p) for x=0,1,2,3,4 assuming p->0.5 and plot the Probability Mass Function.

Use the dbinom function (built into R) which reads as: dbinom(x, size, prob, log = FALSE), where: x = 0,1,2,3,4 (probability of 0:4 seeds germinating) size = 10 (number of seeds planted) prob = 0.5 (given, 50% probability germination/failure) log = FALSE

par(mfrow=c(1,1)) # resets the par function

#Calculate  P(x|p)  for x=0,1,2,3,4 assuming p->0.5  - use the dbinom function
dist <- dbinom(0:4,10,0.5)
dist
## [1] 0.0009765625 0.0097656250 0.0439453125 0.1171875000 0.2050781250
#plot these results
plot(dist, xlab = "n", ylab = "probability", main = "Binomial distribution for x=0,1,2,3,4  assuming p->0.5")

#plot the probability mass function for all of n (size)
dist2 <-dbinom(0:10,10,0.5)
plot(0:10, dist2, type = "b", xlab = "n", ylab = "probability", main = "probability mass for all n (seedling outcomes)")

iii. Plot the Likelihood function given 7 of the 10 seeds germinated. Estimate the mle (maximum likelihood estimate) for p?

#create a binomial distribution for x= 7, n= 10 and p is unknown - so we assign a sequence of 0,1 (failure to germinate, germinate) for all possible options
dist3 <-dbinom(7,10,seq(0,1, by = 0.1))
plot(0:10, dist3, xlab = "n", ylab = "likelihood", main = "Likelihood function for x=7, n=10, and p=unknown")

dist3
##  [1] 0.000000000 0.000008748 0.000786432 0.009001692 0.042467328
##  [6] 0.117187500 0.214990848 0.266827932 0.201326592 0.057395628
## [11] 0.000000000
max(dist3) # the MLE is approx .266 which is the likelihood for 7 [seedlings]
## [1] 0.2668279

Problem 2: Question. Is a coin fair? You sample 18 coin tosses and observe 10 heads. Compare the likelihood function for p in the following two cases: A Bernoulli process that results in c(0,1,1,0,0,1,0,1,0,0,1,1,1,1,0,0,1,1) and a Binomial (k=10,n=18).

i. Compare the pmf for the Bernoulli and Binomial distributions. How are they similar and different?

Bernoulli distribution is the probability distributon over two discrete values of y (0,1) for any fixed value of theta. The binomial is described above in 1.i

Both distributions have the same kernals where one difference is that bernoulli does not have the term \[{{n}\choose{x}}\] while the binomial does (which can be ignored for likelihood).

ii. Plot the likelihood functions in the same figure.

binomDist <- dbinom(10,18,seq(0,1, length.out = 19)) # length.out tells r the number of trails for p (1-19, or 18 trials)
print(binomDist)
##  [1] 0.000000e+00 7.757892e-09 4.891164e-06 1.683042e-04 1.720983e-03
##  [6] 8.859296e-03 2.891444e-02 6.734018e-02 1.194132e-01 1.669235e-01
## [11] 1.865831e-01 1.662890e-01 1.156578e-01 5.988884e-02 2.108204e-02
## [16] 4.207605e-03 3.130345e-04 2.242031e-06 0.000000e+00
binomDistNormal <- binomDist/max(binomDist) # normalize to 1

# create bernoulli likelihood
# this was my first attempt which is INCORRECT, but never figured out quite why
# this is incorrect because it doesn't produce a vector of products:
prod(dbinom(x=(c(0,1,1,0,0,1,0,1,0,0,1,1,1,1,0,0,1,1)),1, seq(0,1, length.out = 19), log=F))
## [1] 0
#this is my second attempt:
#concatenate the results from flips into a variable
bernoulliX <- c(0,1,1,0,0,1,0,1,0,0,1,1,1,1,0,0,1,1) 

# this function is redudant because I couldn't incorporate theta as a sequenced variable, but this works
bernoulliFunc <- function(x) {
  p0 <- prod(dbinom(x,1,0))
  p1 <- prod(dbinom(x,1,0.1))
  p2 <- prod(dbinom(x,1,0.2))
  p3 <- prod(dbinom(x,1,0.3))
  p4 <- prod(dbinom(x,1,0.4))
  p5 <- prod(dbinom(x,1,0.5))
  p6 <- prod(dbinom(x,1,0.6))
  p7 <- prod(dbinom(x,1,0.7))
  p8 <- prod(dbinom(x,1,0.8))
  p9 <- prod(dbinom(x,1,0.9))
  p10 <- prod(dbinom(x,1,1.0))
  return(c(p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10))
}

bernoulliDistNormal <- bernoulliFunc(bernoulliX)/max(bernoulliFunc(bernoulliX)) # normalize to 1

# plot both distributions on same figure
plot(seq(0,1, length.out=11), bernoulliDistNormal, main="Binomial (red) and Bernoulli (dots) Likelihood Functions", xlab="theta", ylab="likelihood")
lines(seq(0,1, length.out=19),binomDistNormal, col="red")

iii. Are the likelihood functions the same? Why or why not?

Yes the functions are basically the same. They differ in how they are contructed but the informaiton/results is the same because the kernals are the same. The only difference is that n is 1 in bernoulli and many in binomial, still everything is proportional:

Binomial \[P(k) = {n\choose k} p^k (1-p)^{(n-k)}\]

Bernoulli \[ p(x) = p^x(1-p)^{(1-x)}\]

iv. Plot confidence limits on the likelihood plot using a likelihood set with alpha =0.1 (see L2.4.2). Is the coin fair?

confidence limit of the maximum to a ratio to a lower value. A very confident confidence limit ratio also makes a very narrow range

Problem 3: Now instead of observing only (k=10,n=18) for a single coin, you are interested in testing the population of coins produced at a US mint. You sample 10 coins, 10 times each. You observe the following number of heads:
6 3 7 3 6 4 4 5 3 8
Plot the likelihood function with the confidence limits as in Problem 2. iv. Is the mint producing fair coins?

noHeads <- (c(6,3,7,3,6,4,4,5,3,8))
# this function is redudant because I couldn't incorporate theta as a sequenced variable, but this works
headsFunc <- function(x) {
  p0 <- prod(dbinom(x,10,0))
  p1 <- prod(dbinom(x,10,0.1))
  p2 <- prod(dbinom(x,10,0.2))
  p3 <- prod(dbinom(x,10,0.3))
  p4 <- prod(dbinom(x,10,0.4))
  p5 <- prod(dbinom(x,10,0.5))
  p6 <- prod(dbinom(x,10,0.6))
  p7 <- prod(dbinom(x,10,0.7))
  p8 <- prod(dbinom(x,10,0.8))
  p9 <- prod(dbinom(x,10,0.9))
  p10 <- prod(dbinom(x,10,1.0))
  return(c(p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10))
}

headsNormalized <- headsFunc(noHeads)/max(headsFunc(noHeads))

plot(seq(0,1, length.out=11), headsNormalized, type = "b", main=" Likelihood Function for 10 coins flipped 10x", xlab="theta", ylab="likelihood")

#errbars(headsNormalized, add=T) # Add error bars confidence intevals

Yes, the likelihood indicates that the coins are concentrated around 0.5, or nearly equal likelihood for heads or tails

Problem 4: Make a figure with 6 panels that show the likelihoods for the following samples (seedlings, seeds): (3,4),(6,8),(12,16),(24,32),(300,400), and (3000,4000).

# as x and n increase, the likelihood lessens and range / spread decreases

par(mfrow=c(2,3))
dist34 <- dbinom(3,4,seq(from=0, to=1, by=0.01))
plot(dist34)
dist68 <- dbinom(6,8,seq(from=0, to=1, by=0.01))
plot(dist68)
dist1216 <- dbinom(12,16,seq(from=0, to=1, by=0.01))
plot(dist1216)
dist2432 <- dbinom(24,32,seq(from=0, to=1, by=0.01))
plot(dist2432)
dist300400 <- dbinom(300,400,seq(from=0, to=1, by=0.01))
plot(dist300400)
dist30004000 <- dbinom(3000,4000,seq(from=0, to=1, by=0.01))
plot(dist30004000)

par(mfrow=c(1,1)) # resets the par function

Problem 5: Let’s simulate seedling counts in forest quadrats. Let’s imagine that we have a field experiment where we are measuring seedlings in forest gaps vs closed canopy.

i. Simulate data that represent 20 quadrats sampled from a forest gap using a lamba (Poisson parameter) of 10 and 20 quadrats from beneath forest canopy with a lamba of 5.

# use rpois to simulate Poisson  distribution: rpois(n, lambda)
quadratsGaps <- rpois(20,10) 
quadratsCanopy <- rpois(20,5)

par(mfrow=c(1,2))
plot(quadratsGaps)
plot(quadratsCanopy)

par(mfrow=c(1,1)) # resets the par function

ii. Plot the likelihood function for the lamba parameters for the gap and forest plots. Estimate the mle for the gap and forest plots.

###likelihood function for a poisson
#e^(-lamda*theta) # taken from page 151 in L2.3

e<-2.71828
thetaPois <- seq(0,1, length.out = 20)
thetaPois
##  [1] 0.00000000 0.05263158 0.10526316 0.15789474 0.21052632 0.26315789
##  [7] 0.31578947 0.36842105 0.42105263 0.47368421 0.52631579 0.57894737
## [13] 0.63157895 0.68421053 0.73684211 0.78947368 0.84210526 0.89473684
## [19] 0.94736842 1.00000000
lambdaGap <- 10
lambdaCanopy <- 5

poisGap <- e^(-lambdaGap*thetaPois)
poisCanopy <- e^(-lambdaCanopy*thetaPois)

par(mfrow=c(1,2))
plot(poisGap, main = "Gap")
plot(poisCanopy, main = "Canopy")

#interprete the MLE from steepness (point of inflection?) tbd
# we can estmimate it by eye, but i future wll need numerical minimizer

par(mfrow=c(1,1)) # resets the par function


####junk
#par(mfrow=c(1,2))
#barplot(table(quadratsGaps), main = "gaps")
#barplot(table(quadratsCanopy), main = "canopy")
#par(mfrow=c(1,1)) # resets the par function

junk code

###**i.     Compare the pmf for the Bernoulli   and Binomial    distributions. How  are they    similar and different?**

###pmf requires p to be known
###assign probability, theta and 1-theta
#thetaHeads <- 10/18
#thetaTails <- 1-thetaHeads

###Create a binomial distribution for k=10, n=18, and p is unknown
#binomDist <-dbinom(1:18,18,thetaHeads) 
#plot(1:18, binomDist, xlab = "# of flips", ylab = "pmf", main = "Binomial distribution for #k=10, n=18, and p inferred from 10 heads in 18 flips")


###Bernouli
###needs work

###same theta
#bernouliDist <-dbinom(1:18,18,thetaHeads) 
#plot(1:18, bernouliDist, xlab = "# of flips", ylab = "pmf")