## ----knitRPreparations,include=FALSE------------------------------------------ library(knitr) opts_chunk$set(tidy=FALSE) ## ----echo = FALSE, results = "hide"------------------------------------------- ## preparation: set option withSweave to TRUE require(distrTEst) require(distrEx) require(distrTeach) require(distrMod) distroptions(withSweave = TRUE) options("newDevice" = TRUE) ## ----exam1, eval = TRUE, fig.width=8.0, fig.height=4.5------------------------ require(distr) N <- Norm(mean = 2, sd = 1.3) P <- Pois(lambda = 1.2) Z <- 2*N + 3 + P Z plot(Z) p(Z)(0.4) q(Z)(0.3) ## in RStudio or Jupyter IRKernel, use q.l(.)(.) instead of q(.)(.) Zs <- r(Z)(50) Zs ## ----DiscrDist, eval = TRUE--------------------------------------------------- D <- DiscreteDistribution(supp = c(1,5,7,21), prob = c(0.1,0.1,0.6,0.2)) D plot(D) ## ----AbscDist, eval = TRUE---------------------------------------------------- AC <- AbscontDistribution(d = function(x) exp(-abs(x)^3), withStand = TRUE) AC plot(AC) ## ----examLis, eval = TRUE----------------------------------------------------- library(distr) M1 <- UnivarMixingDistribution(Norm(), Pois(lambda=1), Norm(), withSimplify = FALSE) M2 <- UnivarMixingDistribution(M1, Norm(), M1, Norm(), withSimplify = FALSE) M2 ## ----warningArith, eval = TRUE------------------------------------------------ A1 <- Norm(); A2 <- Unif() A1 + A2 ## ----examdcP, eval = TRUE----------------------------------------------------- decomposePM(Norm()) decomposePM(Binom(2,0.3)-Binom(5,.4)) decomposePM(UnivarLebDecDistribution(Norm(),Binom(2,0.3)-Binom(5,.4), acWeight = 0.3)) ## ----examflat, eval = TRUE---------------------------------------------------- D1 <- Norm() D2 <- Pois(1) D3 <- Binom(1,.4) D4 <- UnivarMixingDistribution(D1,D2,D3, mixCoeff = c(0.4,0.5,0.1), withSimplify = FALSE) D <- UnivarMixingDistribution(D1,D4,D1,D2, mixCoeff = c(0.4,0.3,0.1,0.2), withSimplify = FALSE) D D0<-flat.mix(D) D0 ## ----arithmetic, eval = TRUE-------------------------------------------------- A1 <- Norm(); A2 <- Unif() d(sin(A1 + A2))(0.1) d(sin(A1 + A2))(0.1) sin(A1 + A2) ## ----arith2v1, eval = TRUE---------------------------------------------------- A1 <- Norm(); A2 <- Unif() A1A2 <- A1*A2 plot(A1A2) ## ----arith2v2, eval = TRUE---------------------------------------------------- A12 <- 1/(A2 + .3) plot(A12) ## ----arith2v3, eval = TRUE---------------------------------------------------- B <- Binom(5,.2)+1 A1B <- A1^B plot(A1B, xlim=c(-3,3)) ## ----arith2v4----------------------------------------------------------------- plot(1.2^A1) ## ----arith2v5----------------------------------------------------------------- plot(B^A1) ## ----Hub, eval = TRUE--------------------------------------------------------- H <- Huberize(Norm(),lower=-1,upper=2) plot(H) ## ----Trun, eval = TRUE-------------------------------------------------------- T <- Truncate(Norm(),lower=-1,upper=2) plot(T) ## ----Min1, eval = TRUE-------------------------------------------------------- M1 <- Maximum(Unif(0,1), Minimum(Unif(0,1), Unif(0,1))) plot(M1) ## ----Min2, eval = TRUE-------------------------------------------------------- M2 <- Minimum(Exp(4),4) plot(M2) ## ----Min3, eval = TRUE-------------------------------------------------------- M3 <- Minimum(Norm(2,2), Pois(3)) plot(M3) ## ----TruncExtr, eval = TRUE--------------------------------------------------- N <- Norm() TN <- Truncate(N, 20,22) r(TN)(20) ## simulation accurate although : p(N)(20, lower.tail = FALSE) ## prob that N>=20 ## ----qrex, eval = TRUE-------------------------------------------------------- B <- Binom(5,0.5) p(B)(3) p.l(B)(3) q(B)(.5) q.r(B)(0.5) ## ----probHN, eval = TRUE------------------------------------------------------ B0 <- as(Binom(5,0.5),"DiscreteDistribution") ## coercion necessary: ## otherwise slot "prob" of B0 will be returned prob(B0) HN <- Huberize(N, -2,1) prob(HN) ## ----makeAC, eval = TRUE------------------------------------------------------ par(mfrow=c(2,3)) plot(makeAbscontDistribution(Nbinom(5,.5)),mfColRow=FALSE) plot(makeAbscontDistribution(HN),mfColRow=FALSE) par(mfrow=c(1,1)) ## ----getLowUp, eval = TRUE---------------------------------------------------- getLow(Nbinom(5,0.5)) getUp(Nbinom(5,0.5)) getLow(Norm(5,0.5)) getUp(Norm(5,0.5)) ## ----cauchy1, eval = TRUE----------------------------------------------------- plot(Cauchy()) ## ----cauchy2, eval = TRUE----------------------------------------------------- plot(Cauchy(),xlim=c(-4,4)) ## ----plotex1, eval = TRUE, fig.width=8.0, fig.height=4.5---------------------- plot(Binom(size = 4, prob = 0.3)) ## ----plotex2, eval = TRUE, fig.width=8.0, fig.height=4.5---------------------- plot(Binom(size = 4, prob = 0.3), do.points = FALSE, verticals = FALSE) ## ----plotex3, eval = TRUE, fig.width=8.0, fig.height=4.5---------------------- plot(Binom(size = 4, prob = 0.3), main = TRUE, inner = FALSE, cex.main = 1.6, tmar = 6) ## ----plotex4, eval = TRUE, fig.width=8.0, fig.height=4.5---------------------- plot(Binom(size = 4, prob = 0.3), cex.points = 1.2, pch = 20, lwd = 2) ## ----plotex5, eval = TRUE, fig.width=8.0, fig.height=4.5---------------------- B <- Binom(size = 4, prob = 0.3) plot(B, col="red", col.points = "green", main = TRUE, col.main="blue", col.sub = "orange", sub = TRUE, cex.sub = 0.6, col.inner = "brown") ## ----plotex6, eval = TRUE, fig.width=8.0, fig.height=4.5---------------------- plot(Nbinom(size = 4,prob = 0.3), cex.points = 1.2, pch.u = 20, pch.a = 10) ## ----plotex7, eval = TRUE, fig.width=8.0, fig.height=4.5---------------------- plot(Chisq(), log = "xy", ngrid = 100) ## ----plotex8, eval = TRUE, fig.width=8.0, fig.height=4.5---------------------- plot(Norm(), lwd=3, col = "red", ngrid = 200, lty = 3, las = 2) ## ----plotex9, eval = TRUE, fig.width=8.0, fig.height=4.5---------------------- plot(Norm(), panel.first = grid(), main = "my Distribution: %A", inner = list(expression(paste(lambda, "-density of %C(%P)")), "CDF", "Pseudo-inverse with param's %N"), sub = "this plot was correctly generated on %D", cex.inner = 0.9, cex.sub = 0.8) ## ----plotex10, eval = TRUE, fig.width=8.0, fig.height=4.5--------------------- Ch <- Chisq(); setgaps(Ch, exactq = 3) plot(Ch, cex = 1.2, pch.u = 20, pch.a = 10, col.points = "green", col.vert = "red") ## ----plotex11, eval = TRUE, fig.width=8.0, fig.height=8.0--------------------- layout(matrix(c(1,3,2,3), nrow=2)) plot(N, mfColRow = FALSE) ## ----plotex12, eval = TRUE, fig.width=8.0, fig.height=8.0--------------------- layout(matrix(c(rep(1,6),2,2,3,3,4,4,5,5,5,6,6,6), nrow=3, byrow=TRUE)) plot(HN, mfColRow = FALSE, to.draw.arg=c("p","d.c","p.c","q.c", "p.d","q.d")) ## ----simulate, eval = TRUE---------------------------------------------------- X <- Simulation() seed(X) <- setRNG() simulate(X) Data(X)[1:10] ## ----expectation, eval = TRUE------------------------------------------------- D4 <- LMCondDistribution(theta = 1) D4 # corresponds to Norm(cond, 1) N <- Norm(mean = 2) E(D4, cond = 1) E(D4, cond = 1, useApply = FALSE) E(as(D4, "UnivariateCondDistribution"), cond = 1) E(D4, function(x){x^2}, cond = 2) E(D4, function(x){x^2}, cond = 2, useApply = FALSE) E(N, function(x){x^2}) E(as(N, "UnivariateDistribution"), function(x){x^2}, useApply = FALSE) # crude Monte-Carlo E(D4, function(x, cond){cond*x^2}, cond = 2, withCond = TRUE) E(D4, function(x, cond){cond*x^2}, cond = 2, withCond = TRUE, useApply = FALSE) E(N, function(x){2*x^2}) E(as(N, "UnivariateDistribution"), function(x){2*x^2}, useApply = FALSE) # crude Monte-Carlo Y <- 5 * Binom(4, .25) - 3 Y E(Y) ## ----expectationlow, eval = TRUE---------------------------------------------- E(Cauchy(), low=3, upp=5) var(Cauchy(), low=3, upp=5) ## ----expectation2, eval = TRUE------------------------------------------------ E(N, function(x)x^2) E(N, function(x)x^2, lowerTruncQuantile = 1e-5) var(Cauchy(), low =3, upperTruncQuantile = 1e-5, IQR.fac = 10) var(Cauchy(), low =3, upperTruncQuantile = 1e-10, IQR.fac = 20) ## ----var, eval = FALSE-------------------------------------------------------- # var <- function(x , ...) # {dots <- list(...) # if(hasArg(y)) y <- dots$"y" # na.rm <- ifelse(hasArg(na.rm), dots$"na.rm", FALSE) # if(!hasArg(use)) # use <- ifelse (na.rm, "complete.obs","all.obs") # else use <- dots$"use" # if(hasArg(y)) # stats::var(x = x, y = y, na.rm = na.rm, use) # else # stats::var(x = x, y = NULL, na.rm = na.rm, use) # } ## ----MCEstimator, eval = TRUE------------------------------------------------- library(distrMod) x <- rgamma(50, scale = 0.5, shape = 3) G <- GammaFamily(scale = 1, shape = 2) negLoglikelihood <- function(x, Distribution){ res <- -sum(log(Distribution@d(x))) names(res) <- "Negative Log-Likelihood" return(res) } MCEstimator(x = x, ParamFamily = G, criterion = negLoglikelihood) ## ----censPoisFamilyDef, eval = FALSE------------------------------------------ # ## search interval for reasonable parameters # startPar <- function(x,...) c(.Machine$double.eps,max(x)) # # ## what to do in case of leaving the parameter domain # makeOKPar <- function(param) {if(param<=0) return(.Machine$double.eps) # return(param)} ## ----PoisFamilyDef, eval = FALSE---------------------------------------------- # setClass("PoisFamily", contains = "L2ParamFamily") ## ----NormLocationFamily, eval = FALSE----------------------------------------- # setClass("NormLocationFamily", contains = "L2LocationFamily") ## ----L2ScaleFamily, eval = FALSE---------------------------------------------- # setMethod("validParameter", signature(object = "L2ScaleFamily"), # function(object, param, tol=.Machine$double.eps){ # if(is(param,"ParamFamParameter")) # param <- main(param) # if(!all(is.finite(param))) return(FALSE) # if(length(param)!=1) return(FALSE) # return(param > tol)}) ## ----GammaFamilyModify, eval = FALSE------------------------------------------ # setMethod("modifyModel", signature(model = "GammaFamily", # param = "ParamFamParameter"), # function(model, param, ...){ # M <- modifyModel(as(model, "L2ParamFamily"), param = param, # .withCall = FALSE) # M@L2derivSymm <- FunSymmList(OddSymmetric(SymmCenter = # prod(main(param))), # NonSymmetric()) # class(M) <- class(model) # return(M) # }) ## ----MLEstimator, eval = TRUE------------------------------------------------- MLEstimator(x = x, ParamFamily = G) MDEstimator(x = x, ParamFamily = G, distance = KolmogorovDist) ## ----NormScaleFam, eval = FALSE----------------------------------------------- # setMethod("mleCalc", signature(x = "numeric", PFam = "NormScaleFamily"), # function(x, PFam, ...){ # n <- length(x) # theta <- sqrt((n-1)/n)*sd(x); mn <- mean(distribution(PFam)) # ll <- -sum(dnorm(x, mean=mn, sd = theta, log=TRUE)) # names(ll) <- "neg.Loglikelihood" # crit.fct <- function(sd) # -sum(dnorm(x, mean=mn, sd = sd, log=TRUE)) # param <- ParamFamParameter(name = "scale parameter", # main = c("sd"=theta)) # if(!hasArg(Infos)) Infos <- NULL # return(meRes(x, theta, ll, param, crit.fct, Infos = Infos)) # }) ## ----CIex, eval = TRUE-------------------------------------------------------- require(distrMod) ## some transformation mtrafo <- function(x){ nms0 <- c("scale","shape") nms <- c("shape","rate") fval0 <- c(x[2], 1/x[1]) names(fval0) <- nms mat0 <- matrix( c(0, -1/x[1]^2, 1, 0), nrow = 2, ncol = 2, dimnames = list(nms,nms0)) list(fval = fval0, mat = mat0)} set.seed(124) x <- rgamma(50, scale = 0.5, shape = 3) ## parametric family of probability measures G <- GammaFamily(scale = 1, shape = 2, trafo = mtrafo) ## MLE res <- MLEstimator(x = x, ParamFamily = G) print(res, digits = 4, show.details="maximal") print(res, digits = 4, show.details="medium") print(res, digits = 4, show.details="minimal") ci <- confint(res) print(ci, digits = 4, show.details="maximal") print(ci, digits = 4, show.details="medium") print(ci, digits = 4, show.details="minimal") ## some profiling par(mfrow=c(2,1)) plot(profile(res)) ## ----NormApprox, eval = TRUE-------------------------------------------------- require(distr) N <- Norm(0,1) U <- Unif(0,1) U2 <- U + U U4 <- U2 + U2 U8 <- U4 + U4 U12 <- U4 + U8 NormApprox <- U12 - 6 x <- seq(-4,4,0.001) opar <- par(no.readonly = TRUE) par(mfrow = c(2,1)) plot(x, d(NormApprox)(x), type = "l", xlab = "", ylab = "Density", main = "Exact and approximated density") lines(x, d(N)(x), col = "red") legend("topleft", legend = c("NormApprox", "Norm(0,1)"), fill = c("black", "red")) plot(x, d(NormApprox)(x) - d(N)(x), type = "l", xlab = "", ylab = "\"black\" - \"red\"", col = "darkgreen", main = "Error") lines(c(-4,4), c(0,0)) par(opar) ## ----ConvolutionNormalDistr, eval = TRUE-------------------------------------- require(distr) ## initialize two normal distributions A <- Norm(mean=1, sd=2) B <- Norm(mean=4, sd=3) ## convolution via addition of moments AB <- A+B ## casting of A,B as absolutely continuous distributions ## that is, ``forget'' that A,B are normal distributions A1 <- as(A, "AbscontDistribution") B1 <- as(B, "AbscontDistribution") ## for higher precision we change the global variable ## "TruncQuantile" from 1e-5 to 1e-8 oldeps <- getdistrOption("TruncQuantile") eps <- 1e-8 distroptions("TruncQuantile" = eps) ## support of A1+B1 for FFT convolution is ## [q(A1)(TruncQuantile), ## q(B1)(TruncQuantile, lower.tail = FALSE)] ## convolution via FFT AB1 <- A1+B1 ############################# ## plots of the results ############################# par(mfrow=c(1,3)) low <- q(AB)(1e-15) upp <- q(AB)(1e-15, lower.tail = FALSE) x <- seq(from = low, to = upp, length = 10000) ## densities plot(x, d(AB)(x), type = "l", lwd = 5) lines(x , d(AB1)(x), col = "orange", lwd = 1) title("Densities") legend("topleft", legend=c("exact", "FFT"), fill=c("black", "orange")) ## cdfs plot(x, p(AB)(x), type = "l", lwd = 5) lines(x , p(AB1)(x), col = "orange", lwd = 1) title("CDFs") legend("topleft", legend=c("exact", "FFT"), fill=c("black", "orange")) ## quantile functions x <- seq(from = eps, to = 1-eps, length = 1000) plot(x, q(AB)(x), type = "l", lwd = 5) lines(x , q(AB1)(x), col = "orange", lwd = 1) title("Quantile functions") legend("topleft", legend=c("exact", "FFT"), fill=c("black", "orange")) ## Since the plots of the results show no ## recognizable differencies, we also compute ## the total variation distance of the densities ## and the Kolmogorov distance of the cdfs ## total variation distance of densities total.var <- function(z, N1, N2){ 0.5*abs(d(N1)(z) - d(N2)(z)) } dv <- integrate(total.var, lower=-Inf, upper=Inf, rel.tol=1e-8, N1=AB, N2=AB1) cat("Total variation distance of densities:\t") print(dv) # 4.25e-07 ### meanwhile realized in package "distrEx" ### as TotalVarDist(N1,N2) ## Kolmogorov distance of cdfs ## the distance is evaluated on a random grid z <- r(Unif(Min=low, Max=upp))(1e5) dk <- max(abs(p(AB)(z)-p(AB1)(z))) cat("Kolmogorov distance of cdfs:\t", dk, "\n") # 2.03e-07 ### meanwhile realized in package "distrEx" ### as KolmogorovDist(N1,N2) ## old distroptions distroptions("TruncQuantile" = oldeps) ## ----ComparisonFFTandRtoDPQ.R, eval = TRUE------------------------------------ require(distr) ################################ ## Comparison 1 - FFT and RtoDPQ ################################ N1 <- Norm(0,3) N2 <- Norm(0,4) rnew1 <- function(n) r(N1)(n) + r(N2)(n) X <- N1 + N2 # exact formula -> N(0,5) Y <- N1 + as(N2, "AbscontDistribution") # appoximated with FFT Z <- new("AbscontDistribution", r = rnew1) # appoximated with RtoDPQ # density-plot x <- seq(-15,15,0.01) plot(x, d(X)(x), type = "l", lwd = 3, xlab = "", ylab = "density", main = "Comparison 1", col = "black") lines(x, d(Y)(x), col = "yellow") lines(x, d(Z)(x), col = "red") legend("topleft", legend = c("Exact", "FFT-Approximation", "RtoDPQ-Approximation"), fill = c("black", "yellow", "red")) ############################################ ## Comparison 2 - "Exact" Formula and RtoDPQ ############################################ B <- Binom(size = 6, prob = 0.5) * 10 N <- Norm() rnew2 <- function(n) r(B)(n) + r(N)(n) Y <- B + N # "exact" formula Z <- new("AbscontDistribution", r = rnew2) # appoximated with RtoDPQ # density-plot x <- seq(-5,65,0.01) plot(x, d(Y)(x), type = "l", xlab = "", ylab = "density", main = "Comparison 2", col = "black") lines(x, d(Z)(x), col = "red") legend("topleft", legend = c("Exact", "RtoDQP-Approximation"), fill = c("black", "red")) ## ----StationaryRegressorDistr, eval = TRUE------------------------------------ require(distr) ## Approximation of the stationary regressor ## distribution of an AR(1) process ## X_t = phi X_{t-1} + V_t ## where V_t i.i.d N(0,1) and phi\in(0,1) ## We obtain ## X_t = \sum_{j=1}^\infty phi^j V_{t-j} ## i.e., X_t \sim N(0,1/(1-phi^2)) phi <- 0.5 ## casting of V as absolutely continuous distributions ## that is, ``forget'' that V is a normal distribution V <- as(Norm(), "AbscontDistribution") ## for higher precision we change the global variable ## "TruncQuantile" from 1e-5 to 1e-8 oldeps <- getdistrOption("TruncQuantile") eps <- 1e-8 distroptions("TruncQuantile" = eps) ## Computation of the approximation ## H=\sum_{j=1}^n phi^j V_{t-j} ## of the stationary regressor distribution ## (via convolution using FFT) H <- V n <- 15 ## may take some time ### switch off warnings [would be issued due to ### very unequal variances...] old.warn <- getOption("warn") options("warn" = -1) for(i in 1:n){Vi <- phi^i*V; H <- H + Vi } options("warn" = old.warn) ## the stationary regressor distribution (exact) X <- Norm(sd=sqrt(1/(1-phi^2))) ############################# ## plots of the results ############################# par(mfrow=c(1,3)) low <- q(X)(1e-15) upp <- q(X)(1e-15, lower.tail = FALSE) x <- seq(from = low, to = upp, length = 10000) ## densities plot(x, d(X)(x),type = "l", lwd = 5) lines(x , d(H)(x), col = "orange", lwd = 1) title("Densities") legend("topleft", legend=c("exact", "FFT"), fill=c("black", "orange")) ## cdfs plot(x, p(X)(x),type = "l", lwd = 5) lines(x , p(H)(x), col = "orange", lwd = 1) title("CDFs") legend("topleft", legend=c("exact", "FFT"), fill=c("black", "orange")) ## quantile functions x <- seq(from = eps, to = 1-eps, length = 1000) plot(x, q(X)(x),type = "l", lwd = 5) lines(x , q(H)(x), col = "orange", lwd = 1) title("Quantile functions") legend( "topleft", legend=c("exact", "FFT"), fill=c("black", "orange")) ## Since the plots of the results show no ## recognizable differencies, we also compute ## the total variation distance of the densities ## and the Kolmogorov distance of the cdfs ## total variation distance of densities total.var <- function(z, N1, N2){ 0.5*abs(d(N1)(z) - d(N2)(z)) } dv <- integrate(f = total.var, lower = -Inf, upper = Inf, rel.tol = 1e-7, N1=X, N2=H) cat("Total variation distance of densities:\t") print(dv) # ~ 5.0e-06 ### meanwhile realized in package "distrEx" ### as TotalVarDist(N1,N2) ## Kolmogorov distance of cdfs ## the distance is evaluated on a random grid z <- r(Unif(Min=low, Max=upp))(1e5) dk <- max(abs(p(X)(z)-p(H)(z))) cat("Kolmogorov distance of cdfs:\t", dk, "\n") # ~2.5e-06 ### meanwhile realized in package "distrEx" ### as KolmogorovDist(N1,N2) ## old distroptions distroptions("TruncQuantile" = oldeps) ## ----destructive, eval = TRUE------------------------------------------------- ########################################################## ## Demo: Instructive destructive example ########################################################## require(distr) ## package "distr" encourages ## consistency but does not ## enforce it---so in general ## d o n o t m o d i f y ## slots d,p,q,r! N <- Norm() B <- Binom() N@d <- B@d plot(N, lwd = 3) ## ----SimulateandEstimate, eval = TRUE, fig.width=8.8, fig.height=3.8---------- require(distrTEst) ### also loads distrSim sim <- new("Simulation", seed = setRNG(), distribution = Norm(mean = 0, sd = 1), filename="sim_01", runs = 1000, samplesize = 30) contsim <- new("Contsimulation", seed = setRNG(), distribution.id = Norm(mean = 0, sd = 1), distribution.c = Norm(mean = 0, sd = 9), rate = 0.1, filename="contsim_01", runs = 1000, samplesize = 30) simulate(sim) simulate(contsim) sim summary(contsim) plot(contsim) ## ----elist, eval = TRUE------------------------------------------------------- require(distrTEst) psim <- function(theta,y,m0){ mean(pmin(pmax(-m0, y - theta), m0)) } mestimator <- function(x, m = 0.7) { uniroot(f = psim, lower = -20, upper = 20, tol = 1e-10, y = x, m0 = m, maxiter = 20)$root } result.id.mean <- evaluate(sim, mean) result.id.mest <- evaluate(sim, mestimator) result.id.median <- evaluate(sim, median) result.cont.mean <- evaluate(contsim, mean) result.cont.mest <- evaluate(contsim, mestimator) result.cont.median <- evaluate(contsim, median) elist <- EvaluationList(result.cont.mean, result.cont.mest, result.cont.median) elist summary(elist) plot(elist, cex = 0.7, las = 2) ## ----mest, eval = TRUE-------------------------------------------------------- result.cont.mest ## ----elist.summary, eval = TRUE----------------------------------------------- summary(elist) ## ----Expectation, eval = TRUE------------------------------------------------- require("distrEx") # Example id <- function(x) x sq <- function(x) x^2 # Expectation and Variance of Binom(6,0.5) B <- Binom(6, 0.5) E(B, id) E(B, sq) - E(B, id)^2 # Expectation and Variance of Norm(1,1) N <- Norm(1, 1) E(N, id) E(N, sq) - E(N, id)^2 ## ----nFoldConvolution, eval = TRUE-------------------------------------------- ########################################################## ## Demo: n-fold convolution of absolutely continuous ## probability distributions ########################################################## require(distr) if(!isGeneric("convpow")) setGeneric("convpow", function(D1,...) standardGeneric("convpow")) ########################################################## ## Function for n-fold convolution ## -- absolute continuous distribution -- ########################################################## ##implentation of Algorithm 3.4. of # Kohl, M., Ruckdeschel, P., Stabla, T. (2005): # General purpose convolution algorithm for distributions # in S4-Classes by means of FFT. # Technical report, Feb. 2005. Also available in # http://www.uni-bayreuth.de/departments/math/org/mathe7/ # /RUCKDESCHEL/pubs/comp.pdf setMethod("convpow", signature(D1 = "AbscontDistribution"), function(D1, N){ if((N < 1)||(!identical(floor(N), N))) stop("N has to be a natural greater than 0") m <- getdistrOption("DefaultNrFFTGridPointsExponent") ##STEP 1 lower <- ifelse((q(D1)(0) > - Inf), q(D1)(0), q(D1)(getdistrOption("TruncQuantile"))) upper <- ifelse((q(D1)(1) < Inf), q(D1)(1), q(D1)(getdistrOption("TruncQuantile"), lower.tail = FALSE)) ##STEP 2 M <- 2^m h <- (upper-lower)/M if(h > 0.01) warning(paste("Grid for approxfun too wide, ", "increase DefaultNrFFTGridPointsExponent", sep="")) x <- seq(from = lower, to = upper, by = h) p1 <- p(D1)(x) ##STEP 3 p1 <- p1[2:(M + 1)] - p1[1:M] ##STEP 4 ## computation of DFT pn <- c(p1, numeric((N-1)*M)) fftpn <- fft(pn) ##STEP 5 ## convolution theorem for DFTs pn <- Re(fft(fftpn^N, inverse = TRUE)) / (N*M) pn <- (abs(pn) >= .Machine$double.eps)*pn i.max <- N*M-(N-2) pn <- c(0,pn[1:i.max]) dn <- pn / h pn <- cumsum(pn) ##STEP 6(density) ## density x <- c(N*lower,seq(from = N*lower+N/2*h, to = N*upper-N/2*h, by=h),N*upper) dnfun1 <- approxfun(x = x, y = dn, yleft = 0, yright = 0) ##STEP 7(density) standardizer <- sum(dn[2:i.max]) + (dn[1]+dn[i.max+1]) / 2 dnfun2 <- function(x) dnfun1(x) / standardizer ##STEP 6(cdf) ## cdf with continuity correction h/2 pnfun1 <- approxfun(x = x+0.5*h, y = pn, yleft = 0, yright = pn[i.max+1]) ##STEP 7(cdf) pnfun2 <- function(x) pnfun1(x) / pn[i.max+1] ## quantile with continuity correction h/2 yleft <- ifelse(((q(D1)(0) == -Inf)| (q(D1)(0) == -Inf)), -Inf, N*lower) yright <- ifelse(((q(D1)(1) == Inf)| (q(D1)(1) == Inf)), Inf, N*upper) w0 <- options("warn") options(warn = -1) qnfun1 <- approxfun(x = pnfun2(x+0.5*h), y = x+0.5*h, yleft = yleft, yright = yright) qnfun2 <- function(x){ ind1 <- (x == 0)*(1:length(x)) ind2 <- (x == 1)*(1:length(x)) y <- qnfun1(x) y <- replace(y, ind1[ind1 != 0], yleft) y <- replace(y, ind2[ind2 != 0], yright) return(y) } options(w0) rnew = function(N) apply(matrix(r(e1)(n*N), ncol=N), 1, sum) return(new("AbscontDistribution", r = rnew, d = dnfun1, p = pnfun2, q = qnfun2)) }) ## initialize a normal distribution A <- Norm(mean=0, sd=1) ## convolution power N <- 10 ## convolution via FFT AN <- convpow(as(A,"AbscontDistribution"), N) ## ... for the normal distribution , 'convpow' has an "exact" ## method by version 1.9 so the as(.,.) is needed to ## see how the algorithm above works ## convolution exact AN1 <- Norm(mean=0, sd=sqrt(N)) ## plots of the results eps <- getdistrOption("TruncQuantile") par(mfrow=c(1,3)) low <- q(AN1)(eps) upp <- q(AN1)(eps, lower.tail = FALSE) x <- seq(from = low, to = upp, length = 10000) ## densities plot(x, d(AN1)(x), type = "l", lwd = 5) lines(x , d(AN)(x), col = "orange", lwd = 1) title("Densities") legend("topleft", legend=c("exact", "FFT"), fill=c("black", "orange")) ## cdfs plot(x, p(AN1)(x), type = "l", lwd = 5) lines(x , p(AN)(x), col = "orange", lwd = 1) title("CDFs") legend("topleft", legend=c("exact", "FFT"), fill=c("black", "orange")) ## quantile functions x <- seq(from = eps, to = 1-eps, length = 1000) plot(x, q(AN1)(x), type = "l", lwd = 5) lines(x , q(AN)(x), col = "orange", lwd = 1) title("Quantile functions") legend("topleft", legend = c("exact", "FFT"), fill = c("black", "orange"))