##### Exercise 3.2: gamma distribution 

loglik.gamma<-function(v,y)
 {
 n<-length(y)
 mu<-mean(y)

 -(-n*lgamma(v)+v*sum(log(y))+n*v*(log(v/mu)-1))

    #### Return minus the log-likelihood,
    #### because optim will minimise this
 }


gamplot<-function(y,max=15000,inc=200)
 {
 # y = data vector
 # max = x axis maximum
 # inc = histogram bin increment


 # Use hist function to get bin frequencies
 carhist<-hist(y,breaks=seq(0,max(y+inc),inc),include.lowest=T,freq=F,plot=F)

 ## Set up axes
  plot(1,0,xlim=c(0,max),ylim=c(0,max(carhist$density)),type="n",
    ylab="f(y)",xlab="Claim size")


 # Draw histogram
  for(i in 1:(length(carhist$breaks)-1))
    if(carhist$breaks[i+1]<max)
        rect(carhist$breaks[i],0,carhist$breaks[i+1],
        carhist$density[i],border=F,col="gray")

 # Compute Gamma MLEs
 v<-optim(1,loglik.gamma,lower=0.01,y=y)$par
 mu <- mean(y)

 # Gamma plot
 x<-seq(0.01,max,inc/10)
 f<-(1/(x*gamma(v)))*(((x*v/mu)^v))*exp(-x*v/mu)  # Gamma density
 lines(x,f)

 }

 gamplot(y=claimcst0[claimcst0>0],max=15000,inc=200)


##### Exercise 3.3: inverse Gaussian distribution 

IGplot<-function(y,max=15000,inc=200)
 {
 # y = data vector
 # max = x axis maximum
 # inc = histogram bin increment


 # Use hist function to get bin frequencies
 carhist<-hist(y,breaks=seq(0,max(y+inc),inc),include.lowest=T,freq=F,plot=F)

 ## Set up axes
  plot(1,0,xlim=c(0,max),ylim=c(0,max(carhist$density)),type="n",
    ylab="f(y)",xlab="Claim size")


 # Draw histogram
  for(i in 1:(length(carhist$breaks)-1))
    if(carhist$breaks[i+1]<max)
        rect(carhist$breaks[i],0,carhist$breaks[i+1],
        carhist$density[i],border=F,col="gray")

 # Compute IG MLEs
 mu<-mean(y)
 sigma2<-(1/length(y))*sum((1/y)-1/mean(y))

 ## IG plot
 x<-seq(0.01,max,inc/10)
 sigma<-sqrt(sigma2)
 f<-(1/(sqrt(2*pi*x^3)*sigma))*exp((-0.5/x)*(((x-mu)/(mu*sigma))^2))
 lines(x,f)

 }

 IGplot(y=claimcst0[claimcst0>0],max=15000,inc=200)