############################################### # CASUALTY ACTUARIAL SOCIETY # # OPEN-SOURCE SOFTWARE COMMITTEE # ############################################### # actuar I Companion Script # # http://opensourcesoftware.casact.org/actuar # ############################################### #Use this script as you read along the article #to quickly execute the examples! #INTRODUCTION library(actuar) #DATA data(dental) data(gdental) dental gdental #PROBABILITY LAWS//FUNCTIONS//PROBABILITY DENSITY FUNCTION dpareto(1.0, shape = 10, scale = 5) #PROBABILITY LAWS//FUNCTIONS//CUMULATIVE DISTRIBUTION FUNCTION ppareto(1.0, shape = 10, scale = 5) #PROBABILITY LAWS//FUNCTIONS//QUANTILE FUNCTION qpareto(.838, shape = 10, scale = 5) #PROBABILITY LAWS//FUNCTIONS//RANDOM DATA r <- rpareto(1000, shape = 10, scale = 5) #PROBABILITY LAWS//FUNCTIONS//RAW MOMENTS m1r <- mpareto(1, shape = 10, scale = 5) m2r <- mpareto(2, shape = 10, scale = 5) variancer <- m2r - m1r^2 stdevr <- variancer^.5 m1r stdevr variancer #PROBABILITY LAWS//FUNCTIONS//LIMITED MOMENTS m1l <- levpareto(1.0, shape = 10, scale = 5, order = 1) m2l <- levpareto(1.0, shape = 10, scale = 5, order = 2) variancel <- m2l - m1l^2 stdevl <- variancel^.5 m1l stdevl variancel rvl <- matrix(c(m1r,m1l,stdevr,stdevl,variancer,variancel),ncol=2,byrow=TRUE) colnames(rvl) <- c("Raw","Limited") rownames(rvl) <- c("Mean","Std Dev","Variance") rvl <- as.table(rvl) rvl #PROBABILITY LAWS//DISTIBUTIONS showgraphs <- function(fun, par, what = c("d", "p", "m", "lev"), xlim) { dist <- switch(fun, trbeta = "TRANSFORMED BETA DISTRIBUTION", genpareto = "GENERALIZED PARETO DISTRIBUTION", burr = "BURR DISTRIBUTION", invburr = "INVERSE BURR DISTRIBUTION", pareto = "PARETO DISTRIBUTION", invpareto = "INVERSE PARETO DISTRIBUTION", llogis = "LOGLOGISTIC DISTRIBUTION", paralogis = "PARALOGISTIC DISTRIBUTION", invparalogis = "INVERSE PARALOGISTIC DISTRIBUTION", trgamma = "TRANSFORMED GAMMA DISTRIBUTION", invtrgamma = "INVERSE TRANSFORMED GAMMA DISTRIBUTION", invgamma = "INVERSE GAMMA DISTRIBUTION", weibull = "WEIBULL DISTRIBUTION", invweibull = "INVERSE WEIBULL DISTRIBUTION", invexp = "INVERSE EXPONENTIAL DISTRIBUTION", pareto1 = "SINGLE PARAMETER PARETO DISTRIBUTION", lgamma = "LOGGAMMA DISTRIBUTION", genbeta = "GENERALIZED BETA DISTRIBUTION", phtype = "PHASE-TYPE DISTRIBUTION", gamma = "GAMMA DISTRIBUTION", exp = "EXPONENTIAL DISTRIBUTION", chisq = "CHI-SQUARE DISTRIBUTION", lnorm = "LOGNORMAL DISTRIBUTION", invgauss = "INVERSE GAUSSIAN DISTRIBUTION", norm = "NORMAL DISTRIBUTION", beta = "BETA DISTRIBUTION", unif = "UNIFORM DISTRIBUTION") if (missing(xlim)) { qf <- match.fun(paste("q", fun, sep = "")) formals(qf)[names(par)] <- par xlim <- c(0, qf(0.999)) } k <- seq.int(4) limit <- seq(0, xlim[2], len = 10) mfrow = c(ceiling(length(what) / 2), 2) op <- par(mfrow = mfrow, oma = c(0, 0, 2, 0)) for (t in what) { f <- match.fun(paste(t, fun, sep = "")) formals(f)[names(par)] <- par main <- switch(t, "d" = "Probability Density Function", "p" = "Cumulative Distribution Function", "m" = "Raw Moments", "lev" = "Limited Expected Value Function", "mgf" = "Moment Generating Function") if (t == "m") plot(k, f(k), type = "l", col = 4, lwd = 2, main = main) else if (t == "lev") plot(limit, f(limit), type = "l", col = 4, lwd = 2, main = main) else if (t == "mgf") curve(f(x), xlim = c(0, 2), col = 4, lwd = 2, main = main) else curve(f(x), xlim = xlim, col = 4, lwd = 2, main = main) title(main = dist, outer = TRUE) } par(op) } showgraphs("invgamma", list(shape = 6, scale = 10)) #PROBABILITY LAWS//DISTIBUTIONS//TRANSFORMED BETA FAMILY showgraphs("trbeta", list(shape1 = 3, shape2 = 4, shape3 = 5, scale = 10)) showgraphs("genpareto", list(shape1 = 10, shape2 = 4, scale = 10)) showgraphs("burr", list(shape1 = 3, shape2 = 4, scale = 10)) showgraphs("invburr", list(shape1 = 3, shape2 = 6, scale = 10)) showgraphs("pareto", list(shape = 10, scale = 10)) showgraphs("invpareto", list(shape = 4, scale = 1), what = c("d", "p")) showgraphs("llogis", list(shape = 6, scale = 10)) showgraphs("paralogis", list(shape = 3, scale = 10)) showgraphs("invparalogis", list(shape = 6, scale = 10)) #PROBABILITY LAWS//DISTRIBUTIONS//TRANSFORMED GAMMA FAMILY showgraphs("trgamma", list(shape1 = 3, shape2 = 1, scale = 10)) showgraphs("invtrgamma", list(shape1 = 3, shape2 = 2, scale = 10)) showgraphs("invgamma", list(shape = 6, scale = 10)) showgraphs("weibull", list(shape = 1.5, scale = 10)) showgraphs("invweibull", list(shape = 6, scale = 10)) showgraphs("invexp", list(rate = 1), what = c("d", "p")) #PROBABILITY LAWS//DISTRIBUTIONS//OTHER showgraphs("pareto1", list(shape = 5, min = 10), xlim = c(0, 50)) showgraphs("lgamma", list(shapelog = 2, ratelog = 5)) showgraphs("genbeta", list(shape1 = 1, shape2 = 2, shape3 = 3, scale = 2)) showgraphs("phtype", list(prob = c(0.5614, 0.4386), rates = matrix(c(-8.64, 0.101, 1.997, -1.095), 2, 2)), what = c("d", "p", "m", "mgf"), xlim = c(0.001, 5)) #PROBABILITY LAWS//DISTRIBUTIONS//INCLUDED IN R showgraphs("gamma", list(shape = 3, rate = 5), what = c("m", "lev", "mgf")) showgraphs("chisq", list(df = 3), what = c("m", "lev", "mgf")) showgraphs("exp", list(rate = 5), what = c("m", "lev", "mgf")) showgraphs("lnorm", list(meanlog = 1, sdlog = 1), what = c("m", "lev")) showgraphs("invgauss", list(nu = 1, lambda = 10), what = c("m", "lev", "mgf"), xlim = c(0, 10)) showgraphs("norm", list(mean = 0, sd = 1), what = c("m", "mgf")) showgraphs("beta", list(shape1 = 1, shape2 = 2), what = c("m", "lev")) showgraphs("unif", list(min = 0, max = 1), what = c("m", "lev", "mgf")) #GROUPED DATA x <- grouped.data(Group = c(0, 25, 50, 100, 150, 250, 500), Line.1 = c(30, 31, 57, 42, 65, 84), Line.2 = c(26, 33, 31, 19, 16, 11)) x #GROUPED DATA//EXTRACTION AND REPLACEMENT x[, 1] x[, -1] x[1:3, ] x[1, 2] <- 22 x x[1, c(2, 3)] <- c(22, 19) x x[1, 1] <- c(0, 20) x x[c(3, 4), 1] <- c(55, 110, 160) x #GROUPED DATA//GROUP MEAN AND HISTOGRAMS mean(x) hist(x[, -3]) #GROUPED DATA//OGIVE CDF <- ogive(x) CDF knots(CDF) plot(CDF) CDF(250) CDF(knots(CDF)) #EMPIRICAL MOMENTS//KTH EMPIRICAL MOMENT emm(dental, order = 1:3) emm(gdental, order = 1:3) #EMPIRICAL MOMENTS//EMPIRICAL LIMITED EXPECTED VALUE lev <- elev(dental) lev glev <- elev(gdental) glev lev(knots(lev)) glev(knots(glev)) plot(lev, type="o", pch = 19) plot(glev, type="o", pch = 19) #OPTIMIZATION//MINIMUM DISTANCE ESTIMATION mde(gdental, pexp, start = list(rate = 1/200), measure = "CvM") mde(gdental, levexp, start = list(rate = 1/200), measure = "LAS") #OPTIMIZATION//MINIMUM DISTANCE ESTIMATION//ISSUES mde(gdental, ppareto, start = list(shape = 3, scale = 600), measure = "CvM") pparetolog <- function(x, logshape, logscale) ppareto(x,exp(logshape),exp(logscale)) pl <- mde(gdental, pparetolog, start = list(logshape = log(3),logscale = log(600)), measure = "CvM") pl exp(pl$estimate) #COVERAGE MODIFICATIONS mcovpdf <- coverage(pdf = dgamma, cdf = pgamma, deductible = 1, limit = 10) mcovcdf <- coverage(cdf = pgamma, deductible = 1, limit = 10) mcovpdf(5, shape = 5, scale = 3) mcovcdf(5, shape = 5, scale = 3) #SOURCES: #C. Dutang, V. Goulet and M. Pigeon (2008). #actuar: An R Package for Actuarial Science. #Journal of Statistical Software, vol. 25, no. 7, 1-37. #URL http://www.jstatsoft.org/v25/i07 #R Development Core Team (2012). #R: A language and environment for statistical computing. #R Foundation for Statistical Computing, Vienna, Austria. #ISBN 3-900051-07-0, URL http://www.R-project.org/.