Nothing
### R code from vignette source 'distrMod.Rnw'
###################################################
### code chunk number 1: library
###################################################
library(distrMod)
distrModoptions(show.details="minimal")
options(prompt = "R> ", continue = "+ ", width = 70,
useFancyQuotes = FALSE)
###################################################
### code chunk number 2: locNorm
###################################################
(N <- NormLocationFamily(mean = 3))
x <- r(N)(20)
###################################################
### code chunk number 3: locNorm1
###################################################
MLEstimator(x,N)
MDEstimator(x,N,distance=CvMDist,
asvar.fct = distrMod:::.CvMMDCovariance)
###################################################
### code chunk number 4: exam1
###################################################
library(distr)
N <- Norm(mean = 2, sd = 1.3)
P <- Pois(lambda = 1.2)
Z <- 2*N + 3 + P
Z
plot(Z, cex.inner = 0.9)
###################################################
### code chunk number 5: exam11
###################################################
p(Z)(0.4)
q(Z)(0.3)
## in RStudio or Jupyter IRKernel, use q.l(.)(.) instead of q(.)(.)
r(Z)(5)
###################################################
### code chunk number 6: expectation
###################################################
library(distrEx)
E(N)
E(P)
E(Z)
E(Z, fun=function(x) sin(x))
###################################################
### code chunk number 7: centralD
###################################################
myD <- AbscontDistribution(d = function(x) exp(-abs(x)^3),
withS = TRUE)
###################################################
### code chunk number 8: locscFam
###################################################
(myFam <- L2LocationScaleFamily(loc = 3, scale = 2,
name = "location and scale family",
centraldistribution = myD))
###################################################
### code chunk number 9: GamFam
###################################################
(G <- GammaFamily(scale = 1, shape = 2))
###################################################
### code chunk number 10: CensPois
###################################################
CensoredPoisFamily <- function(lambda = 1, trunc.pt = 2){
name <- "Censored Poisson family"
distribution <- Truncate(Pois(lambda = lambda),
lower = trunc.pt)
param0 <- lambda
names(param0) <- "lambda"
param <- ParamFamParameter(name = "positive mean",
main = param0,
fixed = c(trunc.pt=trunc.pt))
modifyParam <- function(theta){
Truncate(Pois(lambda = theta),
lower = trunc.pt)}
startPar <- function(x,...) c(.Machine$double.eps,max(x))
makeOKPar <- function(param){
if(param<=0) return(.Machine$double.eps)
return(param)
}
L2deriv.fct <- function(param){
lambda <- main(param)
fct <- function(x){}
body(fct) <- substitute({
x/lambda-ppois(trunc.pt-1,
lambda = lambda,
lower.tail=FALSE)/
ppois(trunc.pt,
lambda = lambda,
lower.tail=FALSE)},
list(lambda = lambda))
return(fct)
}
res <- L2ParamFamily(name = name,
distribution = distribution,
param = param,
modifyParam = modifyParam,
L2deriv.fct = L2deriv.fct,
startPar = startPar,
makeOKPar = makeOKPar)
res@fam.call <- substitute(CensoredPoisFamily(lambda = l,
trunc.pt = t),
list(l = lambda, t = trunc.pt))
return(res)
}
(CP <- CensoredPoisFamily(3,2))
###################################################
### code chunk number 11: plotFam
###################################################
layout(matrix(c(1,1,2,2,3,3,4,4,4,5,5,5), nrow = 2,
byrow = TRUE))
plot(myFam, mfColRow = FALSE, cex.inner = 1,
inner = c("density", "cdf", "quantile function",
"location part", "scale part"))
###################################################
### code chunk number 12: 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)
## })
###################################################
### code chunk number 13: ParaTrafex
###################################################
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,
dimnames = list(nms,nms0))
list(fval = fval0, mat = mat0)
}
(G.MASS <- GammaFamily(scale = 1, shape = 2, trafo = mtrafo))
###################################################
### code chunk number 14: mad-const
###################################################
(mad.const <- 1/q(myD)(0.75))
## in RStudio or Jupyter IRKernel, use q.l(.)(.) instead of q(.)(.)
###################################################
### code chunk number 15: MLE-MASS
###################################################
set.seed(19)
x <- r(distribution(myFam))(50)
mydf <- function(x, loc, scale){
y <- (x-loc)/scale; exp(-abs(y)^3)/scale
}
Med <- median(x)
MAD <- mad(x, constant = mad.const)
c(Med, MAD)
fitdistr(x, mydf, start = list("loc" = Med, "scale" = MAD))
###################################################
### code chunk number 16: MLE-stats4
###################################################
ll <- function(loc,scale){-sum(log(mydf(x,loc,scale)))}
mle(ll, start = list("loc" = Med, "scale" = MAD))
###################################################
### code chunk number 17: MCEstimator
###################################################
negLoglikelihood <- function(x, Distribution){
res <- -sum(log(Distribution@d(x)))
names(res) <- "Negative Log-Likelihood"
return(res)
}
MCEstimator(x = x, ParamFamily = myFam,
criterion = negLoglikelihood)
###################################################
### code chunk number 18: PoisFamilyDef (eval = FALSE)
###################################################
## setClass("PoisFamily", contains = "L2ParamFamily")
###################################################
### code chunk number 19: NormLocationFamily (eval = FALSE)
###################################################
## setClass("NormLocationFamily", contains = "L2LocationFamily")
###################################################
### code chunk number 20: 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)
## })
###################################################
### code chunk number 21: mlecalcEx (eval = FALSE)
###################################################
## setMethod("mleCalc", signature(x = "numeric",
## PFam = "NormLocationScaleFamily"),
## function(x, PFam){
## n <- length(x)
## c(mean(x), sqrt((n-1)/n)*sd(x))
## })
###################################################
### code chunk number 22: MLEstimator
###################################################
MLEstimator(x = x, ParamFamily = myFam)
###################################################
### code chunk number 23: MDEstimatorK
###################################################
MDEstimator(x = x, ParamFamily = myFam,
distance = KolmogorovDist)
###################################################
### code chunk number 24: 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))
## })
###################################################
### code chunk number 25: CIex
###################################################
set.seed(19)
y <- rgamma(50, scale = 3, shape = 2)
(MDest <- MDEstimator(x = y, ParamFamily = G.MASS,
distance = CvMDist))
fitdistr(x = y, densfun = "gamma",
start = list("shape" = estimate(MDest)[1],
"rate" = estimate(MDest)[2]))
(res <- MLEstimator(x = y, ParamFamily = G.MASS))
(ci <- confint(res))
###################################################
### code chunk number 26: profile
###################################################
par(mfrow=c(2,1))
plot(profile(res))
###################################################
### code chunk number 27: cleanup
###################################################
options("prompt"=">")
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.