# Example of maximum likelihood estimation for "n" in a model of # i.i.d. binomial data (with p fixed at 1/2). # The likelihood function. Takes as arguments a vector of data points, # assumed to be i.i.d., and values for n and p. Returns the probability # of the data set if data points have the binomial(n,p) distribution. likelihood <- function (data,n,p) prod(dbinom(data,n,p)) # Plot the likelihood function for a given data set, a vector of n values, # and a given p (default 1/2). plot_likelihood <- function (data,nvals,p=1/2) { lik <- numeric(length(nvals)) for (i in 1:length(nvals)) lik[i] <- likelihood(data,nvals[i],p) plot(nvals,lik,pch=20) } # Estimate n from data by maximum likelihood, with p fixed (default 1/2). estimate_n <- function (data, p=1/2) { n <- max(data)-1 last_likelihood <- 0 while (likelihood(data,n+1,p) >= last_likelihood) { last_likelihood <- likelihood(data,n+1,p) n <- n + 1 } n } # Try out these functions. d <- c(13,9,14) # a small data set plot_likelihood(d,1:30) # plot indicates that n = 24 should be the mle print(estimate_n(d)) # we do indeed get 24 as the mle D <- rep(14,500) # a larger data set print(likelihood(D,28,1/2)) # the likelihood for what ought to be the mle - 0! # print(estimate_n(D)) # goes into an infinite loop if actually done # The log likelihood function. Looking at the log should avoid the problem # of underflow to zero. log_likelihood <- function (data,n,p) sum(dbinom(data,n,p,log=TRUE)) # Estimate n from data by maximum likelihood, with p fixed (default 1/2). # Now done using the log likelihood. estimate_n2 <- function (data, p=1/2) { n <- max(data)-1 last_log_likelihood <- -Inf while (log_likelihood(data,n+1,p) >= last_log_likelihood) { last_log_likelihood <- log_likelihood(data,n+1,p) n <- n + 1 } n } # Try out the version using log likelihood. d <- c(13,9,14) print(estimate_n(d)) # we still get 24 as the mle on this data D <- rep(14,500) print(log_likelihood(D,28,1/2)) # the log likelihood is a good finite value print(estimate_n2(D)) # we find the correct mle of 28