Nothing
inst/doc/distr.R
In distrDoc: Documentation for 'distr' Family of R Packages
## ----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"))
Try the distrDoc package in your browser
Any scripts or data that you put into this service are public.
distrDoc documentation built on Jan. 31, 2024, 3:05 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## ----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"))
Try the distrDoc package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.