1setClass("distr", representation(x = "vector", name = "character", parameters = "numeric", sd = "numeric", n = "numeric", loglik = "numeric"))
2setClass("distrCollection", representation(distr = "list"))
3setMethod("[", signature(x = "distrCollection", i = "ANY"), function(x, i, drop = missing) {
4    x@distr[i]
5})
6setMethod("show", signature(object = "distrCollection"), function(object) {
7    distrList = object@distr
8    cat("\n")
9    for (i in seq(along = distrList)) {
10        temp = distrList[[i]]
11        cat("\n")
12        cat("fitted distribution is", temp@name, ":\n")
13        print(temp@parameters)
14        cat("\n")
15    }
16})
17setMethod("summary", signature(object = "distrCollection"), function(object) {
18    numDist = length(object@distr)
19    gofMatrix = data.frame(matrix(nrow = numDist, ncol = 3))
20    names(gofMatrix) = c("Distribution", "A", "p.value")
21    cat("\n------ Fitted Distribution and estimated parameters ------\n")
22    for (i in seq(along = object@distr)) {
23        distrObj = object@distr[[i]]
24        x = distrObj@x
25        distribution = distrObj@name
26        parameters = distrObj@parameters
27        statistic = NA
28        p.value = NA
29        temp = .myADTest(x, distribution)
30        try(statistic <- as.numeric(temp$statistic), silent = TRUE)
31        try(p.value <- as.numeric(temp$p.value), silent = TRUE)
32        gofMatrix[i, ] = c(distribution, as.numeric(statistic), as.numeric(p.value))
33        cat("\n")
34        cat("fitted distribution is", distribution, ":\n")
35        print(parameters)
36    }
37    cat("\n")
38    cat("\n------ Goodness of Fit - Anderson Darling Test ------\n")
39    cat("\n")
40    gofMatrixPrint = gofMatrix
41    gofMatrixPrint[, 2] = signif(as.numeric(gofMatrixPrint[, 2]), 4)
42    gofMatrixPrint[, 3] = signif(as.numeric(gofMatrixPrint[, 3]), 4)
43    print(gofMatrixPrint)
44})
45distribution = function(x, distribution = "weibull", start, ...) {
46    #if (!require(MASS, quietly = TRUE))
47     #   stop("Package MASS needs to be installed!")
48    if (is.character(distribution))
49        distribution = tolower(distribution)
50    allDistr = c("beta", "cauchy", "chi-squared", "exponential", "f", "gamma", "geometric", "log-normal", "logistic", "negative binomial", "normal", "poisson",
51        "t", "weibull")
52    if (distribution %in% allDistr)
53        distrVec = distribution
54    else distrVec = c("normal")
55    if (identical(distribution, "all"))
56        distrVec = allDistr
57    if (identical(distribution, "quality"))
58        distrVec = c("normal", "log-normal", "exponential", "weibull")
59    distrColl = new("distrCollection")
60    for (i in seq(along = distrVec)) {
61        fit = new("distr")
62        temp = fitdistr(x, distrVec[i], ...)
63        fit@x = x
64        fit@name = distrVec[i]
65        fit@parameters = temp$estimate
66        fit@sd = temp$sd
67        fit@loglik = temp$loglik
68        distrColl@distr[i] = fit
69    }
70    return(distrColl)
71}
72setGeneric("plot", function(x, y, ...) standardGeneric("plot"))
73setMethod("plot", signature(x = "distr"), function(x, y, main, xlab, xlim, ylim, ylab, line.col, line.width, box = TRUE, ...) {
74    object = x
75    xVals = object@x
76    parameters = object@parameters
77    lq = NULL
78    uq = NULL
79    y = NULL
80    if (missing(line.col))
81        line.col = "red"
82    if (missing(line.width))
83        line.width = 1
84    if (missing(main))
85        main = object@name
86    if (missing(xlab))
87        xlab = "x"
88    if (missing(ylab))
89        ylab = "Density"
90    distr = object@name
91    qFun = .charToDistFunc(distr, type = "q")
92    dFun = .charToDistFunc(distr, type = "d")
93    adTestStats = .myADTest(xVals, distr)
94    print(adTestStats)
95    if (class(adTestStats) == "adtest") {
96        A = adTestStats$statistic
97        p = adTestStats$p.value
98    }
99    else {
100        A = NA
101        p = NA
102    }
103    histObj = hist(xVals, plot = FALSE)
104    if (missing(xlim)) {
105        lq = do.call(qFun, c(list(1e-04), as.list(parameters)))
106        uq = do.call(qFun, c(list(0.9999), as.list(parameters)))
107        xlim = range(lq, uq, xVals)
108    }
109    xPoints = seq(xlim[1], xlim[2], length = 200)
110    yPoints = do.call(dFun, c(list(xPoints), as.list(parameters)))
111    if (missing(ylim)) {
112        ylim = range(0, histObj$density, yPoints)
113    }
114    hist(xVals, freq = FALSE, xlab = xlab, ylab = ylab, xlim = xlim, ylim = ylim, main = main, ...)
115    lines(xPoints, yPoints, col = line.col, lwd = line.width)
116    abline(h = 0)
117    legend("topright", c(paste(c(names(parameters), "A", "p"), ": ", c(format(parameters, digits = 3), format(A, digits = 3), format(p, digits = 3))), sep = " "),
118        inset = 0.02)
119    if (box) {
120        box()
121    }
122})
123.xylimits = function(distrCollection, lowerquantile = 0.001, upperquantile = 0.999) {
124    x = NULL
125    y = NULL
126    for (i in seq(along = distrCollection@distr)) {
127        object = distrCollection@distr[[i]]
128        xValues = object@x
129        parameters = object@parameters
130        distr = object@name
131        qFun = .charToDistFunc(distr, type = "q")
132        dFun = .charToDistFunc(distr, type = "d")
133        lq = do.call(qFun, c(list(lowerquantile), as.list(parameters)))
134        uq = do.call(qFun, c(list(upperquantile), as.list(parameters)))
135        x = range(xValues, x, lq, uq)
136        histObj = hist(xValues, plot = FALSE)
137        xPoints = seq(x[1], x[2], length = 200)
138        yPoints = do.call(dFun, c(list(xPoints), as.list(parameters)))
139        y = range(y, 0, histObj$density, yPoints)
140    }
141    invisible(list(xlim = x, ylim = y))
142}
143setGeneric("plot", function(x, y, ...) standardGeneric("plot"))
144setMethod("plot", signature(x = "distrCollection"), function(x, y, xlab, ylab, xlim, ylim, line.col, line.width, ...) {
145    y = NULL
146    object = x
147    distrList = object@distr
148    numDist = length(object@distr)
149    numColWin = ceiling(numDist/2)
150    if (missing(xlim))
151        xlim = .xylimits(object)$xlim
152    if (missing(ylim))
153        ylim = .xylimits(object)$ylim
154    if (missing(line.col))
155        line.col = "red"
156    if (missing(line.width))
157        line.width = 1
158    lapply(distrList, plot, xlim = xlim, ylim = ylim, line.col = line.col, line.width = line.width, ...)
159    cat(paste("Total of", numDist, "plots created"))
160    cat("\n")
161    cat(paste("Use par(mfrow = c(2,", numColWin, ") to see all of them!", sep = ""))
162    cat("\n")
163})
164qqPlot <- function (x, y, confbounds = TRUE, alpha, main, xlab, ylab, xlim, ylim, border = "red", bounds.col = "black", bounds.lty = 1,
165    start, ...)
166{
167    DB = FALSE
168    parList = list(...)
169    if (is.null(parList[["col"]]))
170        parList$col = 1:2
171    if (is.null(parList[["pch"]]))
172        parList$pch = 19
173    if (is.null(parList[["lwd"]]))
174        parList$lwd = 1
175    if (is.null(parList[["cex"]]))
176        parList$cex = 1
177
178    #if (!require(MASS))
179    #    stop("Package MASS needs to be installed!")
180
181    if (class(x) == "distrCollection") {
182        distList = x@distr
183        for (i in 1:length(distList)) {
184            d = distList[[i]]
185            do.call(qqPlot, c(list(x = d@x, y = d@name), parList))
186        }
187        invisible()
188    }
189    if (missing(y))
190        y = "normal"
191	if(missing(alpha))
192		alpha = 0.05
193    if (alpha <=0 || alpha >=1)
194            stop(paste("alpha should be between 0 and 1!"))
195    if (missing(main))
196        main = paste("Q-Q Plot for", deparse(substitute(y)),
197            "distribution")
198    if (missing(xlab))
199        xlab = paste("Quantiles for", deparse(substitute(x)))
200    if (missing(ylab))
201        ylab = paste("Quantiles from", deparse(substitute(y)),
202            "distribution")
203    if (is.numeric(y)) {
204        cat("\ncalling (original) qqplot from namespace stats!\n")
205        return(stats::qqplot(x, y, ...))
206    }
207    qFun = NULL
208    theoretical.quantiles = NULL
209    xs = sort(x)
210    distribution = tolower(y)
211    distWhichNeedParameters = c("weibull", "logistic", "gamma",
212        "exponential", "f", "geometric", "chi-squared", "negative binomial",
213        "poisson")
214
215
216	# new
217	threeParameterDistr = c("weibull3", "lognormal3", "gamma3")
218	threeParameter = distribution %in% threeParameterDistr
219	if(threeParameter) distribution = substr(distribution, 1, nchar(distribution)-1)
220	# end new
221
222	if (is.character(distribution)) {
223        qFun = .charToDistFunc(distribution, type = "q")
224        if (is.null(qFun))
225            stop(paste(deparse(substitute(y)), "distribution could not be found!"))
226    }
227    theoretical.probs = ppoints(xs)
228
229    xq = NULL
230    yq = quantile(xs, prob = c(0.25, 0.75))
231    dots <- list(...)
232    if (TRUE) {
233        if (DB)
234            print("TODO: Pass the estimated parameters correctly")
235        fitList = .lfkp(parList, formals(qFun))
236        fitList$x = xs
237        fitList$densfun = distribution
238        if (!missing(start))
239            fitList$start = start
240        if (DB) {
241            print(fitList)
242            print("Ende")
243        }
244		# new
245		if(!threeParameter){
246        fittedDistr = do.call(fitdistr, fitList)
247        parameter = fittedDistr$estimate
248
249		#save the distribution parameter#
250		thethas = fittedDistr$estimate
251		# save the cariance-covariance matrix
252		varmatrix = fittedDistr$vcov
253		# end of my code
254
255		# new code for three parameter
256		} else {
257			parameter = do.call(paste(".",distribution, "3", sep = ""), list(xs) )    ####
258			threshold = parameter$threshold
259		}
260
261        parameter = .lfkp(as.list(parameter), formals(qFun))
262        params = .lfkp(parList, formals(qFun))
263        parameter = .lfrm(as.list(parameter), params)
264        parameter = c(parameter, params)
265        theoretical.quantiles = do.call(qFun, c(list(c(theoretical.probs)),
266            parameter))
267
268		# new
269		if(!threeParameter){
270		# array containing names of the distributions, for which conf intervals can be computed
271		confIntCapable = c("exponential", "log-normal", "logistic", "normal", "weibull", "gamma", "beta", "cauchy")
272		getConfIntFun = .charToDistFunc(distribution, type = ".confint")
273		# if possible, compute the conf intervals
274		if(confbounds == TRUE){
275			if(distribution %in% confIntCapable){
276				confInt = getConfIntFun(xs, thethas, varmatrix, alpha)
277			}
278		}# end of my code
279		}
280
281        xq <- do.call(qFun, c(list(c(0.25, 0.75)), parameter))
282        if (DB) {
283            print(paste("parameter: ", parameter))
284            print(xq)
285        }
286    }
287    else {
288        params = .lfkp(parList, formals(qFun))
289        params$p = theoretical.probs
290        theoretical.quantiles = do.call(qFun, params)
291        params$p = c(0.25, 0.75)
292        xq = do.call(qFun, params)
293    }
294
295    params = .lfkp(parList, c(formals(plot.default), par()))
296
297	if(!threeParameter){
298		params$y = theoretical.quantiles
299	}  else {
300		params$y = theoretical.quantiles+threshold
301	}
302	params$x = xs
303    params$xlab = xlab
304    params$ylab = ylab
305    params$main = main
306    if (!(is.null(params$col[1]) || is.na(params$col[1])))
307        params$col = params$col[1]
308    if (!missing(xlim))
309        params$xlim = xlim
310    if (!missing(ylim))
311        params$ylim = ylim
312    params$lwd = 1
313    do.call(plot, params)
314    pParams = params
315    pParams = .lfkp(pParams, list(x = 1, y = 1, col = 1, cex = 1))
316    do.call(points, pParams)
317    params = .lfkp(parList, c(formals(abline), par()))
318    params$a = 0
319    params$b = 1
320    params$col = border
321    do.call(abline, params)
322
323	if(!threeParameter){
324	# plot the confInt if available
325	if(confbounds == TRUE){
326	 if(distribution %in% confIntCapable){
327		params = .lfkp(parList, c(formals(lines), par()))
328		params$x = confInt[[3]]
329		params$y = confInt[[1]]
330		params$col = bounds.col
331		params$lty = bounds.lty
332		do.call(lines, params)
333
334		params$x = confInt[[3]]
335		params$y = confInt[[2]]
336		params$col = bounds.col
337		params$lty = bounds.lty
338		do.call(lines, params)
339		}
340	} #end of my function
341	}
342
343    invisible(list(x = theoretical.quantiles, y = xs, int = params$a,
344        slope = params$b))
345}
346
347ppPlot <- function (x, distribution, confbounds = TRUE, alpha, probs, main, xlab, ylab, xlim, ylim,
348    border = "red", bounds.col = "black", bounds.lty = 1, grid = TRUE, box = TRUE, stats = TRUE, start,
349    ...)
350{
351    DB = FALSE
352    conf.level = 0.95
353    conf.lines = TRUE
354    #if (!require(MASS))
355    #    stop("Package MASS needs to be installed!")
356    if (!(is.numeric(x) | (class(x) == "distrCollection")))
357        stop(paste(deparse(substitute(x)), " needs to be numeric or an object of class distrCollection"))
358    parList = list(...)
359    if (is.null(parList[["col"]]))
360        parList$col = c("black", "red", "gray")
361    if (is.null(parList[["pch"]]))
362        parList$pch = 19
363    if (is.null(parList[["lwd"]]))
364        parList$lwd = 1
365    if (is.null(parList[["cex"]]))
366        parList$cex = 1
367    if (DB)
368        print(parList)
369    qFun = NULL
370    xq = NULL
371    yq = NULL
372    x1 = NULL
373	if(missing(alpha))
374		alpha = 0.05
375    if (alpha <=0 || alpha >=1)
376            stop(paste("alpha should be between 0 and 1!"))
377    if (missing(probs))
378        probs = ppoints(11)
379    else if (min(probs) <= 0 || max(probs) >= 1)
380        stop("probs should be values within (0,1)!")
381    probs = round(probs, 2)
382    if (is.numeric(x)) {
383        x1 <- sort(na.omit(x))
384        if (missing(xlim))
385            xlim = c(min(x1) - 0.1 * diff(range(x1)), max(x1) +
386                0.1 * diff(range(x1)))
387    }
388    if (missing(distribution))
389        distribution = "normal"
390    if (missing(ylim))
391        ylim = NULL
392    if (missing(main))
393        main = paste("Probability Plot for", deparse(substitute(distribution)),
394            "distribution")
395    if (missing(xlab))
396        xlab = deparse(substitute(x))
397    if (missing(ylab))
398        ylab = "Probability"
399    if (class(x) == "distrCollection") {
400        distList = x@distr
401        for (i in 1:length(distList)) {
402            d = distList[[i]]
403            do.call(ppPlot, c(list(x = d@x, distribution = d@name),
404                parList))
405        }
406        invisible()
407    }
408    distWhichNeedParameters = c("weibull", "gamma", "logistic",
409        "exponential", "f", "geometric", "chi-squared", "negative binomial",
410        "poisson")
411	# new
412	threeParameterDistr = c("weibull3", "lognormal3", "gamma3")
413	threeParameter = distribution %in% threeParameterDistr
414	if(threeParameter) distribution = substr(distribution, 1, nchar(distribution)-1)
415	# end new
416
417    if (is.character(distribution)) {
418        qFun = .charToDistFunc(distribution, type = "q")
419        pFun = .charToDistFunc(distribution, type = "p")
420        dFun = .charToDistFunc(distribution, type = "d")
421        if (is.null(qFun))
422            stop(paste(deparse(substitute(y)), "distribution could not be found!"))
423    }
424    dots <- list(...)
425    if (TRUE) {
426        if (DB)
427            print("TODO: Pass the estimated parameters correctly")
428        fitList = .lfkp(parList, formals(qFun))
429        fitList$x = x1
430        fitList$densfun = distribution
431        if (!missing(start))
432            fitList$start = start
433
434		if(!threeParameter){
435        fittedDistr = do.call(fitdistr, fitList)
436        parameter = fittedDistr$estimate
437		#save the distribution parameter#
438		thethas = fittedDistr$estimate
439		# save the cariance-covariance matrix
440		varmatrix = fittedDistr$vcov
441		} else {
442			parameter = do.call(paste(".",distribution, "3", sep = ""), list(x1) )    ####
443			print(parameter[3])
444			threshold = parameter$threshold
445		}
446        parameter = .lfkp(as.list(parameter), formals(qFun))
447        params = .lfkp(parList, formals(qFun))
448        parameter = .lfrm(as.list(parameter), params)
449        print(parameter)
450        parameter = c(parameter, params)
451        if (DB) {
452            print(qFun)
453            print(as.list(parameter))
454            print(list(probs))
455        }
456
457
458		# new
459		if(!threeParameter){
460		# array containing names of the distributions, for which conf intervals can be computed
461		confIntCapable = c("exponential", "log-normal", "logistic", "normal", "weibull", "gamma", "beta", "cauchy")
462		getConfIntFun = .charToDistFunc(distribution, type = ".confint")
463		# if possible, compute the conf intervals
464		if(confbounds == TRUE){
465			if(distribution %in% confIntCapable){
466				confInt = getConfIntFun(x1, thethas, varmatrix, alpha)
467			}
468		}# end of my code
469		}
470
471
472		y = do.call(qFun, c(list(ppoints(x1)), as.list(parameter)))
473        yc = do.call(qFun, c(list(ppoints(x1)), as.list(parameter)))
474        cv = do.call(dFun, c(list(yc), as.list(parameter)))
475        print(cv)
476        axisAtY = do.call(qFun, c(list(probs), as.list(parameter)))
477        yq = do.call(qFun, c(list(c(0.25, 0.75)), as.list(parameter)))
478        xq = quantile(x1, probs = c(0.25, 0.75))
479        if (DB) {
480            print(paste("parameter: ", parameter))
481            print(xq)
482        }
483    }
484    else {
485        params = .lfkp(parList, formals(qFun))
486        params$p = ppoints(x1)
487        y = do.call(qFun, params)
488        params$p = probs
489        axisAtY = do.call(qFun, params)
490        params$p = c(0.25, 0.75)
491        yq = do.call(qFun, params)
492        xq = quantile(x1, probs = c(0.25, 0.75))
493    }
494    params = .lfkp(parList, c(formals(plot.default), par()))
495
496	params$x = x1
497	params$y = y
498    params$xlab = xlab
499    params$ylab = ylab
500    params$main = main
501    params$xlim = xlim
502    params$axes = FALSE
503    params$lwd = 1
504    if (!(is.null(params$col[1]) || is.na(params$col[1])))
505        params$col = params$col[1]
506    do.call(plot, params)
507    pParams = params
508    params = .lfkp(parList, list(cex.main = 1, cex.axis = 1,
509        cex.lab = 1))
510    params$side = 1
511    axisAtX = do.call(axis, params)
512    params$side = 2
513    params$at = axisAtY
514    params$labels = probs
515    params$las = 2
516    do.call(axis, params)
517    if (grid) {
518        params = .lfkp(parList, c(formals(abline), list(lwd = 1,
519            col = 1)))
520        params$h = axisAtY
521        params$v = axisAtX
522        if (!(is.null(params$col[3]) || is.na(params$col[3])))
523            params$col = params$col[3]
524        else params$col = 1
525        if (!(is.null(params$lwd[2]) || is.na(params$lwd[2])))
526            params$lwd = params$lwd[2]
527        else params$lwd = 1
528        do.call(abline, params)
529    }
530    pParams = .lfkp(pParams, list(x = 1, y = 1, col = 1, cex = 1))
531    do.call(points, pParams)
532    params = .lfkp(parList, c(formals(abline), par()))
533    if(!threeParameter){
534		params$a = 0
535	}  else {
536		params$a = -threshold
537	}
538    params$b = 1
539    params$col = border
540    do.call(abline, params)
541
542	if(!threeParameter){
543	# plot the confInt if available
544	if(confbounds == TRUE){
545	 if(distribution %in% confIntCapable){
546		params = .lfkp(parList, c(formals(lines), par()))
547		params$x = confInt[[3]]
548		params$y = confInt[[2]]
549		params$col = bounds.col
550		params$lty = bounds.lty
551		do.call(lines, params)
552
553		params$x = confInt[[3]]
554		params$y = confInt[[1]]
555		params$col = bounds.col
556		params$lty = bounds.lty
557		do.call(lines, params)
558		}
559	} #end of my function
560	}
561    if (box)
562        box()
563    invisible(list(x = x, y = y, int = params$a, slope = params$b))
564}
565