# edps 590Bay
# Fall 2019
# C.J.Anderson
#
#  Some examples of Gamma and Inverse Gamma Distributions
#
#  Requires MCMCpack for inverse Gamma
#
x <- seq(from=0, to=15, length.out=1000)

########################### Gamma, fixed b=2.0, vary a #############################
d2.1 <- dgamma(x,1,2)
d2.2<- dgamma(x,2,2)
d2.3<- dgamma(x,3,2)
d2.5<- dgamma(x,5,2)
d2.10<- dgamma(x,10,2)

plot(x,d2.1,type="l", main="Example of Gamma where Scale b=2.0",ylim=c(0,1), col="blue",lwd=2,ylab="Density")
lines(x,d2.2,type="l",lwd=2,col="green")
lines(x,d2.3,type="l",lwd=2,col="purple")
lines(x,d2.5,type="l",lwd=2,col="orange")
lines(x,d2.10,type="l",lwd=2,col="red")

legend("topright",
  legend = c("a=1","a=2","a=3","a=5","a=10"),
  col = c("blue","green","purple","orange","red"),
  text.col = c("blue","green","purple","orange","red"), 
  title="Shape parameter", 
  lty=c(1,1)
)


########################### Gamma, fixed a=2.0, vary b #############################

d2.05<- dgamma(x,2,.5)
d2.075<- dgamma(x,2,.75)
d2.1 <- dgamma(x,2,1)
d2.3<- dgamma(x,2,3)
d2.10<- dgamma(x,2,10)
d2.15<- dgamma(x,2,15)

plot(x,d2.05,type="l", main="Example of Gamma where Shape a=2",ylim=c(0,4), col="blue",lwd=2)
lines(x,d2.075,type="l",lwd=2,col="purple")
lines(x,d2.1,type="l",lwd=2,col="green")
lines(x,d2.3,type="l",lwd=2,col="orange")
lines(x,d2.10,type="l",lwd=2,col="red")
lines(x,d2.15,type="l",lwd=2,col="grey")
legend("topright",
  legend = c("b=0.5","b=0.75","b=1","b=3","b=10","b=15"),
  col = c("blue","purple","green","orange","red","grey"),
  text.col = c("blue","purple","green","orange","red","grey"),
  title="Scale parameter",
  lty=c(1,1)
)


########################### Gamma, fixed a=b #############################

d.05<- dgamma(x,.5,.5)
d.075<- dgamma(x,.75,.75)
d.1 <- dgamma(x,1,1)
d.2<- dgamma(x,2,2)
d.5<- dgamma(x,5,5)
d.10<- dgamma(x,10,10)
d.20<- dgamma(x,20,20)

plot(x,d.05,type="l", main="Example of Gamma where a=b",ylim=c(0,4), col="blue",lwd=2)
lines(x,d.075,type="l",lwd=2,col="purple")
lines(x,d.1,type="l",lwd=2,col="green")
lines(x,d.2,type="l",lwd=2,col="orange")
lines(x,d.5,type="l",lwd=2,col="red")
lines(x,d.10,type="l",lwd=2,col="grey")
lines(x,d.20,type="l",lwd=2,col="cyan")
legend("topright",
  legend = c("b=0.5","a=b=0.75","a=b=1","a=b=2","a=b=5","a=b=10","a=b=20"),
  col = c("blue","purple","green","orange","red","grey","cyan"),
  text.col = c("blue","purple","green","orange","red","grey","cyan"),
  title="Scale parameter",
  lty=c(1,1)
)

########################### Inverse Gamma, fixed b=1.0, vary a #############################
# rinvgamma(n, shape, rate = 1)
# dinvgamma(x, shape, rate = 1)

library(MCMCpack)


d1.1<- dinvgamma(x,1,1)
d2.1<- dinvgamma(x,2,1)
d3.1<- dinvgamma(x,3,1)
d4.1<- dinvgamma(x,4,1)
d5.1<- dinvgamma(x,5,1)

plot(x,d1.1,type="l", main="Inverse Gamma where Scale b=1.0",ylim=c(0,5), col="blue",lwd=2,ylab="Density")
lines(x,d2.1,type="l",lwd=2,col="purple")
lines(x,d3.1,type="l",lwd=2,col="green")
lines(x,d4.1,type="l",lwd=2,col="orange")
lines(x,d5.1,type="l",lwd=2,col="red")
legend("topright",
  legend = c("a=1","a=2","a=3","a=4","a=5"),
  col = c("blue","purple","green","orange","red"),
  text.col = c("blue","purple","green","orange","red"), 
  title="Shape parameter", 
  lty=c(1,1)
)

########################### Inverse Gamma, fixed b=2.0, vary a #############################

d1.2<- dinvgamma(x,1,2)
d2.2<- dinvgamma(x,2,2)
d3.2<- dinvgamma(x,3,2)
d4.2<- dinvgamma(x,4,2)
d5.2<- dinvgamma(x,5,2)

plot(x,d1.2,type="l", main="Inverse Gamma where Scale b=2.0",ylim=c(0,5), col="blue",lwd=2,ylab="Density")
lines(x,d2.2,type="l",lwd=2,col="purple")
lines(x,d3.2,type="l",lwd=2,col="green")
lines(x,d4.2,type="l",lwd=2,col="orange")
lines(x,d5.2,type="l",lwd=2,col="red")
legend("topright",
  legend = c("a=1","a=2","a=3","a=4","a=5"),
  col = c("blue","purple","green","orange","red"),
  text.col = c("blue","purple","green","orange","red"), 
  title="Shape parameter", 
  lty=c(1,1)
)
########################### InverseGamma, fixed a=2.0, vary b #############################

d2.05 <- dinvgamma(x,2,.5)
d2.075<- dinvgamma(x,2,.75)
d2.1  <- dinvgamma(x,2,1)
d2.3 <- dinvgamma(x,2,3)
d2.10<- dinvgamma(x,2,10)
d2.15<- dinvgamma(x,2,15)

plot(x,d2.05,type="l", main="Inverse Gamma where Shape a=2",ylim=c(0,3), col="blue",lwd=2)
lines(x,d2.075,type="l",lwd=2,col="purple")
lines(x,d2.1,type="l",lwd=2,col="green")
lines(x,d2.3,type="l",lwd=2,col="orange")
lines(x,d2.2,type="l",lwd=2,col="red")
lines(x,d2.10,type="l",lwd=2,col="black")
lines(x,d2.075,type="l",lwd=2,col="grey")
legend("topright",
  legend = c("b=0.5","b=0.75","b=1","b=3","b=10","b=15"),
  col = c("blue","purple","green","orange","red","black","grey"),
  text.col = c("blue","purple","green","orange","red","black","grey"),
  title="Scale parameter",
  lty=c(1,1)
)


#########################################################################
# Bigger a and b
#########################################################################
par(mfrow=c(1,1))
big <- seq(from=0, to=40, length.out=1000)
bigg20.30 <- dinvgamma(big,30,30)
bigg10.100 <- dinvgamma(big,10,100)
bigg30.100 <- dinvgamma(big,30,100)
bigg30.10 <- dinvgamma(big,30,10)
bigg30.70 <- dinvgamma(big,30,70)
bigg50.1000 <- dinvgamma(big,50,1000)

plot(big,bigg20.30,type="l", main="Larger values of a and b", col="blue",lwd=2)
lines(big,bigg10.100,col="purple",lwd=2)
lines(big,bigg30.100,col="green",lwd=2)
lines(big,bigg30.10,col="orange",lwd=2)
lines(big,bigg30.70,col="red",lwd=2)
lines(big,bigg50.1000,col="grey",lwd=2)

legend("topright",
  legend = c("a=20, b=20",
             "a=10, b=100",
			 "a=30, b=100",
			 "a=30, b=10",
			 "a=30  b=70",
			 "a=50, b=1000"),
  col = c("blue","purple","green","orange","red","grey"),
  text.col = c("blue","purple","green","orange","red","grey"),
  title="Parameters",
  lty=c(1,1)
)

# To better see the skew in some

plot(big,bigg30.100,type="l", main="Larger values of a and b", col="green",lwd=2)
lines(big,bigg10.100,col="purple",lwd=2)
lines(big,bigg50.1000,col="grey",lwd=2)

legend("topright",
  legend = c("a=30, b=10",
			 "a=30, b=100",
			 "a=50, b=1000"),
  col = c("green","purple","grey"),
  text.col = c("green","purple","grey"),
  title="Parameters",
  lty=c(1,1)
)