---
title: "Geometric likelihood"
author: "Nancy"
date: "07/01/2025"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## 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$.
```{r geomlikelihood and log-lik}
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:



```{r, figures-side, fig.show="hold", out.width="50%"}

## 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

```{r, figures-side2, fig.show="hold", out.width="50%"}

## 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")}

```

```{r, figures-side3, fig.show="hold", out.width="50%"}

## 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")}

```