Question

The daily rainfall in Japan (measured in millimeters) is modelled using a gamma distribution with parameters α = 0.8 and β = 0.4.

Consider the overall rainfall in 365 days, and use mgfs and their properties to prove that this is Ga (292, 0.4).

1. Use the central limit theorem to approximate the probability that the annual rainfall exceeds 800mm (write down the analytical formula and the code used to calculate the cdf value).

2. Compare the approximate value obtained with the exact counterpart (include the R code used to calculate the cdf value).

3. Use R to produce a graphical representation that compares the true and approximate densities (include both the code and plot).

Answer 1

If is the random variable representing the daily rainfall, then the random variable representing the annual rainfall is

. Here . Also

d)According to Central Limit Theorem (CLT), approximately, the distribution of is

The approximate probability that the annual rainfall exceeds 800 mm is

The R commands for the above probability is

> 1-pnorm(800,mean=730,sd=sqrt(1825))
[1] 0.05065079
> 1-pnorm(1.6385761 )
[1] 0.05065079

e) The exact distribution of is (Sum of IID Gamma distribution is Gamma again)

The exact probability that the annual rainfall exceeds 800 mm is

The R commands for the above probability is

> 1-pgamma(800,shape=365*0.8,rate=0.4)
[1] 0.05388408

f) The graphical representation that compares the true and approximate densities is given below. Blue represents the approximate normal and red represents the exact Gamma PDFs. The probability is the area shown shaded.

The R code given below:

plot(1:1)
dev.new()
cord.x1 <- c(800,seq(800,1000,0.1),1000)
cord.y1 <- c(0,dgamma(seq(800,1000,0.1),shape=365*0.8,rate=0.4),0)
cord.x2 <- c(800,seq(800,1000,0.1),1000)
cord.y2 <- c(0,dnorm(seq(800,1000,0.1),mean=730,sd=sqrt(1825)),0)
curve(dgamma(x,shape=365*0.8,rate=0.4), xlim=c(400,1000),ylim = c(0,0.01), main="Gamma & Normal distribution", lwd=2, col = "blue", xlab="x", ylab = "Density")
curve(dnorm(x,mean=730,sd=sqrt(1825)),lwd=2, col="red",add=TRUE)
# Add the shaded area.
polygon(cord.x1,cord.y1,col='skyblue')
polygon(cord.x2,cord.y2,col='violet')