"dgam"<-
function(x, alpha, beta = 1)
{
        dgamma(x/beta, alpha)/beta
}

"pgam"<-
function(x, alpha, beta = 1)
{
        pgamma(x/beta, alpha)
}


"demo.gamma2"<-
function()
{
        par(mfrow=c(2, 2))
        cat("Some examples of the gamma family", fill = TRUE)
        cat("Scale parameter has been chosen so that the distributions", fill
                 = TRUE)
        cat(" have a variance of one", fill = TRUE)
        alpha <- seq(1, 10,  , 15)
        beta <- 1/sqrt(alpha)
        x <- seq(0.01, 6,  , 40)
        z <- matrix(NA, ncol = 15, nrow = 40)
        pz <- matrix(NA, ncol = 15, nrow = 40)
        for(k in 1:15) {
                z[, k] <- dgam(x, alpha = alpha[k], beta = beta[k])
                pz[, k] <- 1-pgam(x, alpha = alpha[k], beta = beta[k])
        }
        matplot(x, z, type = "l", xlab = "x", ylab = "density")
        title("Some examples of the Gamma family\ndensity functions", cex = 
                0.59999999999999998)
        persp(x, alpha, z, xlab = "x", zlab = "Density", ylab = "Alpha", cex = 
                0.5)
        title("Perspective view with shape parameter", cex = 
                0.59999999999999998)    
        #contour(x, alpha, z, nlevels = 12, labex = 0, xlab = "x", ylab = 
#       "shape parameter alpha")
        matplot(x, pz, type = "l", xlab = "x", ylab = "distribution\nfunction")
        title(" Survival distribution functions", cex = 0.59999999999999998)
        image(x, alpha, z, xlab = "x", ylab = "Alpha")
        lines(alpha * beta, alpha)
        lines((alpha - 1) * beta, alpha, lty = 3)
        title("Contour plot of surface\nLines indicate expected values and modes",
                cex = 0.59999999999999998)
}