HW7

1

library(MASS)
data(cats)
fun <- function(mean,variance){
  s <- variance/mean
  a <- mean^2/variance
  c <- c("a"=a,"s"=s)
  return(c)
}
output <- fun(mean = mean(cats$Hwt), variance = var(cats$Hwt))
hist(cats$Hwt, probability = TRUE)
curve(dgamma(x,shape =output[1],scale = output[2]),add = TRUE)

2

fitdistr(cats$Hwt,"gamma")$loglik

## [1] -325.5476

3

gamma.loglike <- function(a,s) {
  logl<- sum(log(dgamma(cats$Hwt,shape = a, scale = s)))
  return(logl)
}
a <- output[1]
s <- output[2]
gamma.loglike(a,s)

## [1] -325.6886

4

gamma.loglike.vec <- Vectorize(gamma.loglike, vectorize.args = c("a","s"))
curve(gamma.loglike.vec(a = x, s= output[2]), from = 1, to=40)

curve(gamma.loglike.vec(a = output[1], s= x), from = 0.01, to=1)

The peak value appears around -325 to -326.

5

surface <- function(expr,from.x=0,to.x=1,from.y=0,to.y=1,n.x=101,
  n.y=101,...) {
  x.seq <- seq(from=from.x,to=to.x,length.out=n.x)
  y.seq <- seq(from=from.y,to=to.y,length.out=n.y)
  unevaluated.expression <- substitute(expr)
  z <- function(x,y) {
    return(eval(unevaluated.expression,envir=list(x=x,y=y)))
  }
  z.values <- outer(X=x.seq,Y=y.seq,FUN=z)
  z.matrix <- matrix(z.values,nrow=n.x)
  contour(x=x.seq,y=y.seq,z=z.matrix,...)
  invisible(list(x=x.seq,y=y.seq,z=z.matrix, func=z))
}
surface(gamma.loglike.vec(a=x,s=y),from.x=1,to.x=40,from.y=0.01,to.y=1,nlevels=1000)

#points(gamma.loglike(output[1], output[2]))

6

Even though they have the approximately same values, the plots from #4 are 2D graphs and #5 is a 3D one, so they are not compatible.

7

surface <- function(expr,from.x=0,to.x=1,from.y=0,to.y=1,n.x=101,
  n.y=101,...) {
  x.seq <- seq(from=from.x,to=to.x,length.out=n.x)
  y.seq <- seq(from=from.y,to=to.y,length.out=n.y)
  unevaluated.expression <- substitute(expr)
  z <- function(x,y) {
    return(eval(unevaluated.expression,envir=list(x=x,y=y)))
  }
  z.values <- outer(X=x.seq,Y=y.seq,FUN=z)
  z.matrix <- matrix(z.values,nrow=n.x)
  contour(x=x.seq,y=y.seq,z=z.matrix,...)
  invisible(list(x=x.seq,y=y.seq,z=z.matrix, func=z))
}
surface(gamma.loglike.vec(a=x,s=y),from.x=19.4,to.x=20.5,from.y=0.5355,to.y=0.54,nlevels=1000)

8

From the plot, we can see that maximum is convergent to -325.568. Compared with -326, it is almost 0.5 greater. If I keep zomming in, the maxmium may increase, but so far for this plot, the maximum is -325.568.

9

step1 <- function(a,s,x) {
  n <- nrow(x)
  logl <- (a-1)*sum(x) - n*gamma(a) - n*a*log(s) - (1/s)*sum(x)
  return(-logl)
}