Geometric likelihood

Setup

The density is \[ f(x;\theta) = \theta (1-\theta)^{(x-1)}, \quad x=1, 2, \dots, \] although in R it is defined for \(x=0, 1, \dots\). The likelihood and log-likelihod functions for a sample of size \(n\) are \[L(\theta;\underline{x}) = \theta^n (1-\theta)^{s-n}, \quad \ell(\theta;\underline{x}) = n\log\theta +(s-n)\log(1-\theta), \quad s=\Sigma x_i.\] You can check that \(\hat\theta = n/s = 1/\bar x\).

geomlik <- function(theta,s, n){
  theta^n*(1-theta)^(s-n)}

geomllik <- function(theta, s, n){
  log(geomlik(theta,s,n))-max(log(geomlik(theta,s,n)))}

Plots

The plots I showed in class had a fairly weird shape (not looking very parabolic). This depends a lot on the range for the y-axis. Here are a few variations:

## set the parameters, draw a single plot
n <- 10; prob <- 0.5
x <- rgeom(n, prob) + 1 #R definition different from mine
s <- sum(x)

thvals <- seq(0.01,0.99,length=100)

plot(thvals,geomllik(thvals, s, n), type="l", lwd=2, ylim=c(-4,0), ylab="log-likelihood", xlab = expression(theta), main = paste("n = ", n, "true theta =", prob))
#then add some more
for(i in 1:15){
  x <- rgeom(n,prob)+1
  s <- sum(x)
  lines(thvals,geomllik(thvals,s,n), col="gray")}

## change a parameter, add another plot
n <- 10; prob <- 0.2
x <- rgeom(n, prob) + 1 #R definition different from mine
s <- sum(x)

thvals <- seq(0.01,0.99,length=100)

plot(thvals,geomllik(thvals, s, n), type="l", lwd=2, ylim=c(-4,0), ylab="log-likelihood", xlab = expression(theta), main = paste("n = ", n, "true theta =", prob))
#then add some more
for(i in 1:15){
  x <- rgeom(n,prob)+1
  s <- sum(x)
  lines(thvals,geomllik(thvals,s,n), col="gray")}

If you remove “ylim=c(-4,0)” from the plot parameters, you get pictures similar to those I showed in class. You can try editing my R Markdown file and knitting it again.

increasing sample size

## set the parameters, draw a single plot
n <- 100; prob <- 0.5
x <- rgeom(n, prob) + 1 #R definition different from mine
s <- sum(x)

thvals <- seq(0.01,0.99,length=100)

plot(thvals,geomllik(thvals, s, n), type="l", lwd=2, ylim=c(-4,0), ylab="log-likelihood", xlab = expression(theta), main = paste("n = ", n, "true theta =", prob))
#then add some more
for(i in 1:15){
  x <- rgeom(n,prob)+1
  s <- sum(x)
  lines(thvals,geomllik(thvals,s,n), col="gray")}

## change a parameter, add another plot
n <- 100; prob <- 0.2
x <- rgeom(n, prob) + 1 #R definition different from mine
s <- sum(x)

thvals <- seq(0.01,0.99,length=100)

plot(thvals,geomllik(thvals, s, n), type="l", lwd=2, ylim=c(-4,0), ylab="log-likelihood", xlab = expression(theta), main = paste("n = ", n, "true theta =", prob))
#then add some more
for(i in 1:15){
  x <- rgeom(n,prob)+1
  s <- sum(x)
  lines(thvals,geomllik(thvals,s,n), col="gray")}

## set the parameters, draw a single plot
n <- 200; prob <- 0.5
x <- rgeom(n, prob) + 1 #R definition different from mine
s <- sum(x)

thvals <- seq(0.01,0.99,length=100)

plot(thvals,geomllik(thvals, s, n), type="l", lwd=2, ylim=c(-4,0), ylab="log-likelihood", xlab = expression(theta), main = paste("n = ", n, "true theta =", prob))
#then add some more
for(i in 1:15){
  x <- rgeom(n,prob)+1
  s <- sum(x)
  lines(thvals,geomllik(thvals,s,n), col="gray")}

## change a parameter, add another plot
n <- 200; prob <- 0.2
x <- rgeom(n, prob) + 1 #R definition different from mine
s <- sum(x)

thvals <- seq(0.01,0.99,length=100)

plot(thvals,geomllik(thvals, s, n), type="l", lwd=2, ylim=c(-4,0), ylab="log-likelihood", xlab = expression(theta), main = paste("n = ", n, "true theta =", prob))
#then add some more
for(i in 1:15){
  x <- rgeom(n,prob)+1
  s <- sum(x)
  lines(thvals,geomllik(thvals,s,n), col="gray")}