Nothing
R/distdata.R
In R4VN: Health Data Analysis and Publication-Ready Reporting
Defines functions
.r4vn_viewer_distdata
.r4vn_dist_svg
plot.r4vn_distdata
distdata
.r4vn_dist_grid
.r4vn_dist_display_table
.r4vn_dist_quantile_table
.r4vn_dist_range_table
.r4vn_dist_probability_table
.r4vn_dist_support_text
.r4vn_dist_param_text
.r4vn_dist_args
.r4vn_dist_validate
.r4vn_dist_find
.r4vn_dist_registry
.r4vn_dist_safe_var
.r4vn_dist_spec
.r4vn_dist_param
Documented in
distdata
plot.r4vn_distdata
# ============================================================================
# R4VN probability distributions: command engine and teaching calculator
# ============================================================================
.r4vn_dist_param <- function(name, label, default, min = NULL, max = NULL, step = NULL, description = "") {
list(name = name, label = label, default = default, min = min, max = max, step = step, description = description)
}
.r4vn_dist_spec <- function(id, label, type, category, aliases, params,
d, p, q, r, support, moments, history, use, notes = NULL,
formula = NULL) {
list(id = id, label = label, type = type, category = category,
aliases = unique(c(id, aliases)), params = params,
d = d, p = p, q = q, r = r, support = support, moments = moments,
history = history, use = use, notes = notes, formula = formula)
}
.r4vn_dist_safe_var <- function(x) if (is.finite(x) && x >= 0) x else NA_real_
.r4vn_dist_registry <- function() {
P <- .r4vn_dist_param
S <- .r4vn_dist_spec
reg <- list()
reg$normal <- S(
"normal", "Normal (Gaussian)", "continuous", "Continuous - symmetric",
c("norm", "gaussian"),
list(P("mean", "Mean (mu)", 0, -10, 10, .1, "Center of the distribution."),
P("sd", "Standard deviation (sigma)", 1, .05, 10, .05, "Controls spread; must be positive.")),
function(x,a) stats::dnorm(x,a$mean,a$sd), function(x,a) stats::pnorm(x,a$mean,a$sd),
function(pr,a) stats::qnorm(pr,a$mean,a$sd), function(n,a) stats::rnorm(n,a$mean,a$sd),
function(a) c(-Inf,Inf), function(a) c(mean=a$mean,variance=a$sd^2,sd=a$sd),
"The bell-shaped law emerged from work by Abraham de Moivre on the binomial approximation in the 18th century and was developed further by Laplace and Gauss in probability and error theory.",
"Common for approximately symmetric continuous measurements, sampling distributions, measurement error, and many regression assumptions. It should not be chosen only because a variable is continuous.",
"About 68%, 95%, and 99.7% of values lie within 1, 2, and 3 standard deviations of the mean when the distribution is exactly normal.",
"f(x) = exp[-(x-mu)^2/(2 sigma^2)] / (sigma sqrt(2 pi))"
)
reg$lognormal <- S(
"lognormal", "Log-normal", "continuous", "Continuous - positive",
c("lnorm", "log-normal"),
list(P("meanlog", "Mean on log scale", 0, -5, 5, .1, "Mean of log(X)."),
P("sdlog", "SD on log scale", 1, .05, 3, .05, "Standard deviation of log(X).")),
function(x,a) stats::dlnorm(x,a$meanlog,a$sdlog), function(x,a) stats::plnorm(x,a$meanlog,a$sdlog),
function(pr,a) stats::qlnorm(pr,a$meanlog,a$sdlog), function(n,a) stats::rlnorm(n,a$meanlog,a$sdlog),
function(a)c(0,Inf), function(a){m<-exp(a$meanlog+a$sdlog^2/2);v<-(exp(a$sdlog^2)-1)*exp(2*a$meanlog+a$sdlog^2);c(mean=m,variance=v,sd=sqrt(v))},
"The log-normal family became important through 19th- and early-20th-century work on multiplicative variation; it describes a variable whose logarithm is normally distributed.",
"Useful for positive right-skewed quantities produced by multiplicative processes, such as some biological measurements, costs, incubation times, and environmental concentrations.",
"The parameters describe log(X), not X. The median on the original scale is exp(meanlog).",
"log(X) ~ Normal(meanlog, sdlog^2)"
)
reg$uniform <- S(
"uniform", "Uniform", "continuous", "Continuous - bounded",
c("unif"),
list(P("min", "Minimum (a)", 0, -10, 10, .1, "Lower support boundary."), P("max", "Maximum (b)", 1, -10, 20, .1, "Upper support boundary; must exceed a.")),
function(x,a) stats::dunif(x,a$min,a$max), function(x,a) stats::punif(x,a$min,a$max), function(pr,a) stats::qunif(pr,a$min,a$max), function(n,a) stats::runif(n,a$min,a$max),
function(a)c(a$min,a$max), function(a)c(mean=(a$min+a$max)/2,variance=(a$max-a$min)^2/12,sd=(a$max-a$min)/sqrt(12)),
"Uniform models are among the oldest probability models and formalize equal density over a bounded interval.",
"Useful as a simple bounded model, for random-number generation, simulation, and as a non-informative building block. It is usually unrealistic for biological measurements unless equal density is scientifically plausible.",
"Every interval of the same width inside [a,b] has the same probability.", "f(x)=1/(b-a), a <= x <= b"
)
reg$exponential <- S(
"exponential", "Exponential", "continuous", "Continuous - positive",
c("exp"), list(P("rate", "Rate (lambda)", 1, .01, 10, .05, "Event rate per unit time; mean equals 1/rate.")),
function(x,a) stats::dexp(x,a$rate), function(x,a) stats::pexp(x,a$rate), function(pr,a) stats::qexp(pr,a$rate), function(n,a) stats::rexp(n,a$rate),
function(a)c(0,Inf), function(a)c(mean=1/a$rate,variance=1/a$rate^2,sd=1/a$rate),
"The exponential distribution arose naturally from Poisson-process theory and waiting-time problems.",
"Used for positive waiting times when the instantaneous hazard is constant. In survival analysis, that constant-hazard assumption should be considered explicitly.",
"It is memoryless: conditional remaining time has the same distribution regardless of elapsed time.", "f(x)=lambda exp(-lambda x), x >= 0"
)
reg$gamma <- S(
"gamma", "Gamma", "continuous", "Continuous - positive",
c("gam"), list(P("shape", "Shape (alpha)", 2, .05, 15, .05, "Controls skewness and shape."), P("rate", "Rate (beta)", 1, .01, 10, .05, "Inverse scale; mean = shape/rate.")),
function(x,a) stats::dgamma(x,a$shape,rate=a$rate), function(x,a) stats::pgamma(x,a$shape,rate=a$rate), function(pr,a) stats::qgamma(pr,a$shape,rate=a$rate), function(n,a) stats::rgamma(n,a$shape,rate=a$rate),
function(a)c(0,Inf), function(a)c(mean=a$shape/a$rate,variance=a$shape/a$rate^2,sd=sqrt(a$shape)/a$rate),
"The family is tied to Euler's gamma function and became central in waiting-time and Pearson-family models.",
"Useful for positive right-skewed continuous outcomes, waiting times to multiple events, costs, and as a flexible GLM distribution with a log link.",
"R4VN uses shape-rate parameterization. `scale` can be supplied to distdata() as an alternative, where rate = 1/scale.", "f(x)=beta^alpha x^(alpha-1) exp(-beta x) / Gamma(alpha)"
)
reg$beta <- S(
"beta", "Beta", "continuous", "Continuous - bounded",
c("beta"), list(P("shape1", "Shape 1 (alpha)", 2, .05, 15, .05, "Controls behavior near 0 and overall shape."), P("shape2", "Shape 2 (beta)", 2, .05, 15, .05, "Controls behavior near 1 and overall shape.")),
function(x,a) stats::dbeta(x,a$shape1,a$shape2), function(x,a) stats::pbeta(x,a$shape1,a$shape2), function(pr,a) stats::qbeta(pr,a$shape1,a$shape2), function(n,a) stats::rbeta(n,a$shape1,a$shape2),
function(a)c(0,1), function(a){m<-a$shape1/(a$shape1+a$shape2);v<-a$shape1*a$shape2/((a$shape1+a$shape2)^2*(a$shape1+a$shape2+1));c(mean=m,variance=v,sd=sqrt(v))},
"The beta family is based on Euler's beta function and later became a core bounded distribution in probability and Bayesian statistics.",
"Useful for proportions or probabilities that vary continuously between 0 and 1 and as a conjugate prior for a Bernoulli/binomial probability.",
"alpha=beta gives symmetry; values below 1 can create U-shaped densities, while values above 1 often create unimodal densities.", "f(x)=x^(alpha-1)(1-x)^(beta-1)/B(alpha,beta)"
)
reg$chisq <- S(
"chisq", "Chi-square", "continuous", "Continuous - sampling distributions",
c("chi-square", "chi square", "chi2"), list(P("df", "Degrees of freedom", 4, .1, 50, .1, "Controls shape and spread."), P("ncp", "Noncentrality", 0, 0, 25, .1, "0 gives the central chi-square distribution.")),
function(x,a) stats::dchisq(x,a$df,ncp=a$ncp), function(x,a) stats::pchisq(x,a$df,ncp=a$ncp), function(pr,a) stats::qchisq(pr,a$df,ncp=a$ncp), function(n,a) stats::rchisq(n,a$df,ncp=a$ncp),
function(a)c(0,Inf), function(a){if(a$ncp==0)c(mean=a$df,variance=2*a$df,sd=sqrt(2*a$df)) else c(mean=a$df+a$ncp,variance=2*(a$df+2*a$ncp),sd=sqrt(2*(a$df+2*a$ncp)))},
"Karl Pearson introduced the chi-square goodness-of-fit framework around 1900; the distribution also arises as a sum of squared standard normal variables.",
"Used in categorical-data tests, variance inference under normality, likelihood theory, and many asymptotic test statistics.",
"For the central distribution, the mean is df and variance is 2*df. A noncentrality parameter is useful in power calculations.", "X = sum(Z_i^2) for df independent standard normal variables"
)
reg$t <- S(
"t", "Student's t", "continuous", "Continuous - sampling distributions",
c("student", "student-t", "student t"), list(P("df", "Degrees of freedom", 5, .5, 100, .5, "Smaller df gives heavier tails."), P("ncp", "Noncentrality", 0, -10, 10, .1, "0 gives the central t distribution.")),
function(x,a) stats::dt(x,a$df,ncp=a$ncp), function(x,a) stats::pt(x,a$df,ncp=a$ncp), function(pr,a) stats::qt(pr,a$df,ncp=a$ncp), function(n,a) stats::rt(n,a$df,ncp=a$ncp),
function(a)c(-Inf,Inf), function(a){if(a$ncp!=0)return(c(mean=NA_real_,variance=NA_real_,sd=NA_real_));m<-if(a$df>1)0 else NA_real_;v<-if(a$df>2)a$df/(a$df-2) else if(a$df>1) Inf else NA_real_;c(mean=m,variance=v,sd=sqrt(v))},
"William Sealy Gosset published the distribution under the pseudonym 'Student' in 1908 while working at Guinness; Fisher later developed its role in classical inference.",
"Used for inference on means and regression coefficients when a standard error is estimated, especially in small samples under normal-model assumptions.",
"As df increases, the central t distribution approaches the standard normal; small df produce heavier tails.", NULL
)
reg$f <- S(
"f", "F", "continuous", "Continuous - sampling distributions",
c("f-distribution", "f distribution"), list(P("df1", "Numerator df", 5, .5, 50, .5, "Numerator degrees of freedom."), P("df2", "Denominator df", 20, .5, 100, .5, "Denominator degrees of freedom."), P("ncp", "Noncentrality", 0, 0, 25, .1, "0 gives the central F distribution.")),
function(x,a) stats::df(x,a$df1,a$df2,ncp=a$ncp), function(x,a) stats::pf(x,a$df1,a$df2,ncp=a$ncp), function(pr,a) stats::qf(pr,a$df1,a$df2,ncp=a$ncp), function(n,a) stats::rf(n,a$df1,a$df2,ncp=a$ncp),
function(a)c(0,Inf), function(a){if(a$ncp!=0)return(c(mean=NA_real_,variance=NA_real_,sd=NA_real_));m<-if(a$df2>2)a$df2/(a$df2-2) else NA_real_;v<-if(a$df2>4)2*a$df2^2*(a$df1+a$df2-2)/(a$df1*(a$df2-2)^2*(a$df2-4)) else NA_real_;c(mean=m,variance=v,sd=sqrt(v))},
"The F distribution is associated with Ronald A. Fisher's variance-ratio methods and was named F by George W. Snedecor.",
"Used in ANOVA, nested-model tests, variance-ratio problems, and regression/model comparison. Noncentral F distributions are common in power calculations.",
"It is a ratio of two appropriately scaled independent chi-square variables in the central case.", NULL
)
reg$weibull <- S(
"weibull", "Weibull", "continuous", "Continuous - positive",
c("weib"), list(P("shape", "Shape", 1.5, .1, 8, .05, "Determines how hazard changes over time."), P("scale", "Scale", 1, .05, 10, .05, "Sets the time/measurement scale.")),
function(x,a) stats::dweibull(x,a$shape,a$scale), function(x,a) stats::pweibull(x,a$shape,a$scale), function(pr,a) stats::qweibull(pr,a$shape,a$scale), function(n,a) stats::rweibull(n,a$shape,a$scale),
function(a)c(0,Inf), function(a){m<-a$scale*gamma(1+1/a$shape);v<-a$scale^2*(gamma(1+2/a$shape)-gamma(1+1/a$shape)^2);c(mean=m,variance=v,sd=sqrt(v))},
"The distribution is named for Waloddi Weibull, who demonstrated its broad applicability in reliability work in the 20th century.",
"Widely used for survival and reliability data. Shape <1 gives decreasing hazard, shape=1 gives the exponential model, and shape >1 gives increasing hazard.",
"Its flexible hazard pattern makes it useful when constant hazard is implausible but monotonic hazard is reasonable.", NULL
)
reg$logistic <- S(
"logistic", "Logistic", "continuous", "Continuous - symmetric",
c("logis"), list(P("location", "Location", 0, -10, 10, .1, "Center and median."), P("scale", "Scale", 1, .05, 10, .05, "Positive spread parameter.")),
function(x,a) stats::dlogis(x,a$location,a$scale), function(x,a) stats::plogis(x,a$location,a$scale), function(pr,a) stats::qlogis(pr,a$location,a$scale), function(n,a) stats::rlogis(n,a$location,a$scale),
function(a)c(-Inf,Inf), function(a)c(mean=a$location,variance=pi^2*a$scale^2/3,sd=pi*a$scale/sqrt(3)),
"The logistic curve is associated with Pierre-Francois Verhulst's 19th-century population-growth work; the corresponding probability distribution later became important in statistical modeling.",
"Useful as a symmetric heavy-tailed alternative to the normal and fundamental to logistic regression through its cumulative distribution/logit relationship.",
"Compared with a normal distribution of similar spread, the logistic has heavier tails.", NULL
)
reg$cauchy <- S(
"cauchy", "Cauchy", "continuous", "Continuous - symmetric",
c("lorentz"), list(P("location", "Location", 0, -10, 10, .1, "Median and center of symmetry."), P("scale", "Scale", 1, .05, 10, .05, "Positive scale parameter.")),
function(x,a) stats::dcauchy(x,a$location,a$scale), function(x,a) stats::pcauchy(x,a$location,a$scale), function(pr,a) stats::qcauchy(pr,a$location,a$scale), function(n,a) stats::rcauchy(n,a$location,a$scale),
function(a)c(-Inf,Inf), function(a)c(mean=NA_real_,variance=NA_real_,sd=NA_real_),
"Named for Augustin-Louis Cauchy, this heavy-tailed law also appears as the ratio of two independent standard normal variables.",
"Useful as a heavy-tailed model and in robust/Bayesian work, but not as a routine substitute for normality. Its extreme tails materially affect summaries and inference.",
"The ordinary mean and variance do not exist, even though the distribution is symmetric around its location parameter.", NULL
)
reg$laplace <- S(
"laplace", "Laplace (double exponential)", "continuous", "Continuous - symmetric",
c("double exponential"), list(P("location", "Location", 0, -10, 10, .1, "Center and median."), P("scale", "Scale", 1, .05, 10, .05, "Mean absolute deviation scale; must be positive.")),
function(x,a) exp(-abs(x-a$location)/a$scale)/(2*a$scale),
function(x,a) ifelse(x<a$location,.5*exp((x-a$location)/a$scale),1-.5*exp(-(x-a$location)/a$scale)),
function(pr,a) a$location + a$scale*ifelse(pr<.5,log(2*pr),-log(2*(1-pr))),
function(n,a){u<-stats::runif(n);a$location+a$scale*ifelse(u<.5,log(2*u),-log(2*(1-u)))},
function(a)c(-Inf,Inf), function(a)c(mean=a$location,variance=2*a$scale^2,sd=sqrt(2)*a$scale),
"The distribution is associated with Pierre-Simon Laplace and is also called the double-exponential distribution.",
"Useful for symmetric data with a sharp center and heavier tails than the normal; it also underlies least-absolute-deviation ideas and the lasso's common Bayesian prior interpretation.",
"The scale parameter is not the standard deviation; SD = sqrt(2)*scale.", NULL
)
reg$gumbel <- S(
"gumbel", "Gumbel (extreme value type I)", "continuous", "Continuous - extreme values",
c("ev1", "extreme value"), list(P("location", "Location", 0, -10, 10, .1, "Shifts the distribution."), P("scale", "Scale", 1, .05, 10, .05, "Positive spread parameter.")),
function(x,a){z<-(x-a$location)/a$scale;exp(-(z+exp(-z)))/a$scale},
function(x,a){z<-(x-a$location)/a$scale;exp(-exp(-z))},
function(pr,a)a$location-a$scale*log(-log(pr)), function(n,a){u<-stats::runif(n);a$location-a$scale*log(-log(u))},
function(a)c(-Inf,Inf), function(a){eg<-0.5772156649015329;m<-a$location+eg*a$scale;v<-pi^2*a$scale^2/6;c(mean=m,variance=v,sd=sqrt(v))},
"Named for Emil Gumbel, a major developer of extreme-value statistics in the 20th century.",
"Used to model maxima or minima (after an appropriate sign convention), such as extreme environmental measurements or maximum loads.",
"This implementation uses the distribution for maxima with CDF exp(-exp(-(x-location)/scale)).", NULL
)
reg$pareto <- S(
"pareto", "Pareto type I", "continuous", "Continuous - heavy tailed",
c("pareto1"), list(P("shape", "Shape (alpha)", 3, .1, 10, .05, "Tail index; larger values give lighter tails."), P("scale", "Minimum (x_m)", 1, .05, 10, .05, "Lower support boundary.")),
function(x,a){out<-rep(0,length(x));ok<-is.finite(x)&x>=a$scale;out[ok]<-a$shape*a$scale^a$shape/x[ok]^(a$shape+1);out[is.infinite(x)&x>0]<-0;out},
function(x,a){out<-rep(0,length(x));ok<-is.finite(x)&x>=a$scale;out[ok]<-1-(a$scale/x[ok])^a$shape;out[is.infinite(x)&x>0]<-1;out},
function(pr,a)a$scale/(1-pr)^(1/a$shape), function(n,a){u<-stats::runif(n);a$scale/(1-u)^(1/a$shape)},
function(a)c(a$scale,Inf), function(a){m<-if(a$shape>1)a$shape*a$scale/(a$shape-1) else Inf;v<-if(a$shape>2)a$shape*a$scale^2/((a$shape-1)^2*(a$shape-2)) else Inf;c(mean=m,variance=v,sd=sqrt(v))},
"Vilfredo Pareto used power-law ideas in studying wealth and income; the Pareto family later became a standard heavy-tail model.",
"Useful for strongly right-skewed heavy-tailed quantities and threshold exceedances in some contexts. Tail assumptions should be checked because moments may not exist.",
"The mean is infinite when shape <=1 and the variance is infinite when shape <=2.", NULL
)
reg$rayleigh <- S(
"rayleigh", "Rayleigh", "continuous", "Continuous - positive",
c("rayleigh"), list(P("scale", "Scale (sigma)", 1, .05, 10, .05, "Positive scale parameter.")),
function(x,a) ifelse(x<0,0,x/a$scale^2*exp(-x^2/(2*a$scale^2))),
function(x,a) ifelse(x<0,0,1-exp(-x^2/(2*a$scale^2))),
function(pr,a)a$scale*sqrt(-2*log(1-pr)), function(n,a)a$scale*sqrt(-2*log(stats::runif(n))),
function(a)c(0,Inf), function(a){m<-a$scale*sqrt(pi/2);v<-(2-pi/2)*a$scale^2;c(mean=m,variance=v,sd=sqrt(v))},
"Named for Lord Rayleigh; the distribution arises as the magnitude of a two-dimensional vector with independent zero-mean normal components of equal variance.",
"Used for magnitudes, signal amplitudes, wind-speed-like models, and radial error when two orthogonal components are Gaussian.",
"It is positive and right-skewed; the single scale parameter controls both center and spread.", NULL
)
reg$triangular <- S(
"triangular", "Triangular", "continuous", "Continuous - bounded",
c("triangle"), list(P("min", "Minimum (a)", 0, -10, 10, .1, "Lower bound."), P("mode", "Mode (c)", .5, -10, 20, .1, "Most likely value; must lie between min and max."), P("max", "Maximum (b)", 1, -10, 20, .1, "Upper bound.")),
function(x,a){out<-rep(0,length(x));i<-x>=a$min&x<=a$mode;out[i]<-2*(x[i]-a$min)/((a$max-a$min)*(a$mode-a$min));j<-x>a$mode&x<=a$max;out[j]<-2*(a$max-x[j])/((a$max-a$min)*(a$max-a$mode));out},
function(x,a){out<-ifelse(x<a$min,0,ifelse(x>=a$max,1,NA_real_));i<-x>=a$min&x<=a$mode;out[i]<-(x[i]-a$min)^2/((a$max-a$min)*(a$mode-a$min));j<-x>a$mode&x<a$max;out[j]<-1-(a$max-x[j])^2/((a$max-a$min)*(a$max-a$mode));out},
function(pr,a){fc<-(a$mode-a$min)/(a$max-a$min);ifelse(pr<fc,a$min+sqrt(pr*(a$max-a$min)*(a$mode-a$min)),a$max-sqrt((1-pr)*(a$max-a$min)*(a$max-a$mode)))},
function(n,a){u<-stats::runif(n);fc<-(a$mode-a$min)/(a$max-a$min);ifelse(u<fc,a$min+sqrt(u*(a$max-a$min)*(a$mode-a$min)),a$max-sqrt((1-u)*(a$max-a$min)*(a$max-a$mode)))},
function(a)c(a$min,a$max), function(a){m<-(a$min+a$mode+a$max)/3;v<-(a$min^2+a$mode^2+a$max^2-a$min*a$mode-a$min*a$max-a$mode*a$max)/18;c(mean=m,variance=v,sd=sqrt(v))},
"The triangular distribution is a simple bounded model long used in simulation and project-risk settings when only a minimum, mode, and maximum are elicited.",
"Useful for transparent teaching and rough uncertainty modeling when information is limited to plausible bounds and a most likely value; not a substitute for empirical distribution fitting.",
"The mode must lie between the lower and upper bounds.", NULL
)
reg$kumaraswamy <- S(
"kumaraswamy", "Kumaraswamy", "continuous", "Continuous - bounded",
c("kuma"), list(P("shape1", "Shape a", 2, .05, 15, .05, "First shape parameter."), P("shape2", "Shape b", 2, .05, 15, .05, "Second shape parameter.")),
function(x,a) ifelse(x<0|x>1,0,a$shape1*a$shape2*x^(a$shape1-1)*(1-x^a$shape1)^(a$shape2-1)),
function(x,a) ifelse(x<=0,0,ifelse(x>=1,1,1-(1-x^a$shape1)^a$shape2)),
function(pr,a)(1-(1-pr)^(1/a$shape2))^(1/a$shape1), function(n,a){u<-stats::runif(n);(1-(1-u)^(1/a$shape2))^(1/a$shape1)},
function(a)c(0,1), function(a){m<-a$shape2*beta(1+1/a$shape1,a$shape2);e2<-a$shape2*beta(1+2/a$shape1,a$shape2);v<-e2-m^2;c(mean=m,variance=v,sd=sqrt(max(0,v)))},
"P. Kumaraswamy proposed this bounded two-shape distribution in 1980 for hydrologic applications.",
"Useful as a flexible beta-like model on (0,1), particularly in simulation because its CDF and quantile function have simple closed forms.",
"It resembles the beta family in flexibility but uses different shape parameters and is not interchangeable with beta regression parameterizations.", NULL
)
reg$loglogistic <- S(
"loglogistic", "Log-logistic", "continuous", "Continuous - positive",
c("log-logistic", "fisk"), list(P("shape", "Shape", 2.5, .1, 10, .05, "Tail/hazard shape."), P("scale", "Scale", 1, .05, 10, .05, "Median and positive scale.")),
function(x,a){out<-rep(0,length(x));ok<-is.finite(x)&x>0;z<-x[ok]/a$scale;out[ok]<-(a$shape/a$scale)*z^(a$shape-1)/(1+z^a$shape)^2;out[is.infinite(x)&x>0]<-0;out},
function(x,a){out<-rep(0,length(x));ok<-is.finite(x)&x>0;z<-x[ok]/a$scale;out[ok]<-1/(1+z^(-a$shape));out[is.infinite(x)&x>0]<-1;out},
function(pr,a)a$scale*(pr/(1-pr))^(1/a$shape), function(n,a){u<-stats::runif(n);a$scale*(u/(1-u))^(1/a$shape)},
function(a)c(0,Inf), function(a){m<-if(a$shape>1)a$scale*(pi/a$shape)/sin(pi/a$shape) else Inf;e2<-if(a$shape>2)a$scale^2*(2*pi/a$shape)/sin(2*pi/a$shape) else Inf;v<-if(is.finite(e2))e2-m^2 else Inf;c(mean=m,variance=v,sd=sqrt(v))},
"The log-logistic distribution results when log(X) follows a logistic distribution and is also known as the Fisk distribution in some fields.",
"Used for positive duration/survival data when a non-monotone hazard may be plausible, and for some income or size distributions.",
"The scale parameter is the median. The mean exists only for shape >1 and variance only for shape >2.", NULL
)
reg$halfnormal <- S(
"halfnormal", "Half-normal", "continuous", "Continuous - positive",
c("half-normal"), list(P("sd", "Underlying SD (sigma)", 1, .05, 10, .05, "SD of the zero-mean normal variable before taking absolute value.")),
function(x,a) ifelse(x<0,0,2*stats::dnorm(x,0,a$sd)), function(x,a) ifelse(x<0,0,2*stats::pnorm(x/a$sd)-1),
function(pr,a)a$sd*stats::qnorm((pr+1)/2), function(n,a)abs(stats::rnorm(n,0,a$sd)), function(a)c(0,Inf),
function(a){m<-a$sd*sqrt(2/pi);v<-a$sd^2*(1-2/pi);c(mean=m,variance=v,sd=sqrt(v))},
"The half-normal is the distribution of the absolute value of a zero-mean normal variable.",
"Useful for non-negative magnitudes, random-effect scale priors, measurement-error magnitudes, and teaching transformations of normal variables.",
"Its parameter is the SD of the underlying normal, not the SD of the half-normal itself.", NULL
)
# Discrete distributions ---------------------------------------------------
reg$bernoulli <- S(
"bernoulli", "Bernoulli", "discrete", "Discrete - counts/events",
c("bern"), list(P("p", "Probability of success", .5, 0, 1, .01, "Probability that X=1.")),
function(x,a) stats::dbinom(x,1,a$p), function(x,a) stats::pbinom(x,1,a$p), function(pr,a) stats::qbinom(pr,1,a$p), function(n,a) stats::rbinom(n,1,a$p),
function(a)c(0,1), function(a)c(mean=a$p,variance=a$p*(1-a$p),sd=sqrt(a$p*(1-a$p))),
"Named for Jacob Bernoulli, whose work on repeated binary trials was published posthumously in Ars Conjectandi in 1713.",
"Models one binary trial: event/no event, success/failure, diseased/not diseased. Repeated independent Bernoulli trials lead to the binomial distribution.",
"X takes only 0 and 1; p is both P(X=1) and E(X).", "P(X=x)=p^x(1-p)^(1-x), x in {0,1}"
)
reg$binomial <- S(
"binomial", "Binomial", "discrete", "Discrete - counts/events",
c("binom"), list(P("n", "Number of trials (n)", 10, 1, 200, 1, "Number of independent Bernoulli trials."), P("p", "Probability of success", .5, 0, 1, .01, "Success probability for each trial.")),
function(x,a) stats::dbinom(x,a$n,a$p), function(x,a) stats::pbinom(x,a$n,a$p), function(pr,a) stats::qbinom(pr,a$n,a$p), function(nn,a) stats::rbinom(nn,a$n,a$p),
function(a)c(0,a$n), function(a)c(mean=a$n*a$p,variance=a$n*a$p*(1-a$p),sd=sqrt(a$n*a$p*(1-a$p))),
"The binomial law grew from early probability work on repeated trials, especially Jacob Bernoulli's law of large numbers and binomial-trial reasoning.",
"Used for the number of events in a fixed number of independent trials when each trial has the same event probability. Examples include positive tests among n specimens or responders among n independent participants under a simple model.",
"Key assumptions are fixed n, two outcomes per trial, constant p, and independence. Overdispersion or clustering can invalidate the simple binomial model.", "P(X=x)=choose(n,x) p^x (1-p)^(n-x)"
)
reg$poisson <- S(
"poisson", "Poisson", "discrete", "Discrete - counts/events",
c("pois"), list(P("lambda", "Mean/rate (lambda)", 3, .01, 30, .05, "Expected count in the specified exposure interval.")),
function(x,a) stats::dpois(x,a$lambda), function(x,a) stats::ppois(x,a$lambda), function(pr,a) stats::qpois(pr,a$lambda), function(n,a) stats::rpois(n,a$lambda),
function(a)c(0,Inf), function(a)c(mean=a$lambda,variance=a$lambda,sd=sqrt(a$lambda)),
"Simeon Denis Poisson described the distribution in 1837 as a limiting law for rare events.",
"Used for counts over time, person-time, area, or another exposure when events occur independently at a stable rate. Poisson regression extends this idea to rate modeling with covariates and offsets.",
"The basic Poisson distribution has mean equal to variance; extra-Poisson variation suggests overdispersion, clustering, heterogeneity, or another count model.", "P(X=x)=exp(-lambda) lambda^x/x!"
)
reg$geometric <- S(
"geometric", "Geometric", "discrete", "Discrete - waiting counts",
c("geom"), list(P("p", "Success probability", .25, .001, 1, .01, "Probability of success on each independent trial.")),
function(x,a) stats::dgeom(x,a$p), function(x,a) stats::pgeom(x,a$p), function(pr,a) stats::qgeom(pr,a$p), function(n,a) stats::rgeom(n,a$p),
function(a)c(0,Inf), function(a)c(mean=(1-a$p)/a$p,variance=(1-a$p)/a$p^2,sd=sqrt(1-a$p)/a$p),
"The geometric distribution is one of the classical waiting-time distributions arising from Bernoulli trials.",
"Used for the number of failures before the first success under independent trials with constant success probability.",
"R4VN follows R's convention: X counts failures before the first success, so support starts at 0. Some textbooks count the trial number instead and therefore start at 1.", NULL
)
reg$nbinom <- S(
"nbinom", "Negative binomial", "discrete", "Discrete - counts/events",
c("negative binomial", "negbin"), list(P("size", "Size/dispersion", 2, .05, 30, .05, "Positive shape/dispersion parameter."), P("p", "Success probability", .4, .001, .999, .01, "R probability parameter; mean=size*(1-p)/p.")),
function(x,a) stats::dnbinom(x,size=a$size,prob=a$p), function(x,a) stats::pnbinom(x,size=a$size,prob=a$p), function(pr,a) stats::qnbinom(pr,size=a$size,prob=a$p), function(n,a) stats::rnbinom(n,size=a$size,prob=a$p),
function(a)c(0,Inf), function(a){m<-a$size*(1-a$p)/a$p;v<-a$size*(1-a$p)/a$p^2;c(mean=m,variance=v,sd=sqrt(v))},
"The negative-binomial family developed from classical repeated-trial problems and later became important as a Poisson-gamma mixture for overdispersed counts.",
"Used for overdispersed count data and heterogeneous event rates. Negative-binomial regression is a common alternative when Poisson variance is too small.",
"distdata() also accepts `mu` instead of `p`; it converts using p=size/(size+mu).", NULL
)
reg$hypergeom <- S(
"hypergeom", "Hypergeometric", "discrete", "Discrete - sampling without replacement",
c("hypergeometric", "hyper"), list(P("success", "Successes in population", 10, 0, 100, 1, "Number of population items classified as success."), P("failure", "Failures in population", 20, 0, 100, 1, "Number of population items classified as failure."), P("draws", "Sample size", 5, 0, 100, 1, "Number drawn without replacement.")),
function(x,a) stats::dhyper(x,a$success,a$failure,a$draws), function(x,a) stats::phyper(x,a$success,a$failure,a$draws), function(pr,a) stats::qhyper(pr,a$success,a$failure,a$draws), function(n,a) stats::rhyper(n,a$success,a$failure,a$draws),
function(a)c(max(0,a$draws-a$failure),min(a$draws,a$success)), function(a){N<-a$success+a$failure;m<-a$draws*a$success/N;v<-if(N>1)a$draws*(a$success/N)*(a$failure/N)*(N-a$draws)/(N-1) else 0;c(mean=m,variance=v,sd=sqrt(v))},
"The hypergeometric distribution is a classical finite-population model closely connected with combinatorial probability.",
"Used when sampling without replacement from a finite population containing a known number of successes and failures; it also underlies Fisher's exact test conditioning arguments.",
"Unlike the binomial distribution, draws are dependent because sampled items are not replaced.", NULL
)
reg$discreteuniform <- S(
"discreteuniform", "Discrete uniform", "discrete", "Discrete - bounded",
c("discrete uniform", "dunifint"), list(P("min", "Minimum integer", 1, -20, 20, 1, "Smallest possible integer."), P("max", "Maximum integer", 6, -20, 50, 1, "Largest possible integer.")),
function(x,a) ifelse(x==floor(x)&x>=a$min&x<=a$max,1/(a$max-a$min+1),0),
function(x,a) pmax(0,pmin(1,(floor(x)-a$min+1)/(a$max-a$min+1))),
function(pr,a) {
k <- a$max - a$min + 1
ifelse(pr <= 0, a$min, ifelse(pr >= 1, a$max,
pmax(a$min, pmin(a$max, a$min + ceiling(pr * k) - 1))))
},
function(n,a) sample(seq.int(a$min,a$max),n,replace=TRUE), function(a)c(a$min,a$max),
function(a){m<-(a$min+a$max)/2;k<-a$max-a$min+1;v<-(k^2-1)/12;c(mean=m,variance=v,sd=sqrt(v))},
"The discrete-uniform model formalizes equally likely values from a finite set, such as a fair die.",
"Useful for randomization examples, equally likely integer outcomes, and introductory probability teaching.",
"Every integer from min through max has exactly the same probability.", NULL
)
reg$betabinom <- S(
"betabinom", "Beta-binomial", "discrete", "Discrete - overdispersed binomial",
c("beta-binomial", "beta binomial"), list(P("n", "Number of trials", 10, 1, 200, 1, "Number of Bernoulli trials."), P("shape1", "Alpha", 2, .05, 20, .05, "First beta shape parameter."), P("shape2", "Beta", 2, .05, 20, .05, "Second beta shape parameter.")),
function(x,a){ok<-x==floor(x)&x>=0&x<=a$n;out<-rep(0,length(x));out[ok]<-exp(lchoose(a$n,x[ok])+lbeta(x[ok]+a$shape1,a$n-x[ok]+a$shape2)-lbeta(a$shape1,a$shape2));out},
function(x,a){vapply(x,function(xx){k<-0:min(a$n,floor(xx));if(xx<0)0 else if(xx>=a$n)1 else sum(exp(lchoose(a$n,k)+lbeta(k+a$shape1,a$n-k+a$shape2)-lbeta(a$shape1,a$shape2)))},numeric(1))},
function(pr,a){vapply(pr,function(pp){
if (pp <= 0) return(0)
if (pp >= 1) return(a$n)
k<-0:a$n
cdf<-cumsum(exp(lchoose(a$n,k)+lbeta(k+a$shape1,a$n-k+a$shape2)-lbeta(a$shape1,a$shape2)))
hit <- which(cdf >= pp)[1L]
if (is.na(hit)) a$n else k[hit]
},numeric(1))},
function(nn,a){pp<-stats::rbeta(nn,a$shape1,a$shape2);stats::rbinom(nn,a$n,pp)}, function(a)c(0,a$n),
function(a){p<-a$shape1/(a$shape1+a$shape2);rho<-1/(a$shape1+a$shape2+1);m<-a$n*p;v<-a$n*p*(1-p)*(1+(a$n-1)*rho);c(mean=m,variance=v,sd=sqrt(v))},
"The beta-binomial combines a binomial count with beta-distributed heterogeneity in the event probability.",
"Useful for overdispersed binomial counts, clustered binary outcomes, and teaching how heterogeneity in p inflates variance beyond the simple binomial model.",
"The induced intraclass-correlation/overdispersion parameter is rho = 1/(alpha+beta+1).", NULL
)
reg$zip <- S(
"zip", "Zero-inflated Poisson", "discrete", "Discrete - zero-inflated counts",
c("zero-inflated poisson", "zero inflated poisson"), list(P("lambda", "Poisson mean", 3, .01, 30, .05, "Mean of the Poisson count component."), P("zi", "Structural-zero probability", .2, 0, .95, .01, "Probability of an extra structural zero.")),
function(x,a){base<-(1-a$zi)*stats::dpois(x,a$lambda);base[x==0]<-a$zi+(1-a$zi)*exp(-a$lambda);base},
function(x,a) ifelse(x<0,0,a$zi+(1-a$zi)*stats::ppois(x,a$lambda)),
function(pr,a){m0<-a$zi+(1-a$zi)*exp(-a$lambda);vapply(pr,function(pp){if(pp<=m0)0 else stats::qpois((pp-a$zi)/(1-a$zi),a$lambda)},numeric(1))},
function(n,a){z<-stats::rbinom(n,1,a$zi);ifelse(z==1,0,stats::rpois(n,a$lambda))}, function(a)c(0,Inf),
function(a){m<-(1-a$zi)*a$lambda;v<-(1-a$zi)*a$lambda*(1+a$zi*a$lambda);c(mean=m,variance=v,sd=sqrt(v))},
"Zero-inflated count models were developed to represent data generated by a mixture of a structural-zero state and an ordinary count process.",
"Useful when there is a defensible two-process mechanism producing more zeros than a Poisson model can explain, for example structural non-risk plus a count process among those at risk.",
"A large number of zeros alone does not prove zero inflation; negative-binomial or hurdle models may be more appropriate depending on the data-generating process.", NULL
)
reg$zipf <- S(
"zipf", "Zipf (finite)", "discrete", "Discrete - heavy tailed",
c("zipfian"), list(P("exponent", "Exponent (s)", 1.5, .1, 5, .05, "Controls the power-law decay."), P("max", "Maximum rank", 50, 2, 500, 1, "Finite upper rank used by this teaching implementation.")),
function(x,a){k<-seq_len(a$max);den<-sum(k^(-a$exponent));ifelse(x==floor(x)&x>=1&x<=a$max,x^(-a$exponent)/den,0)},
function(x,a){k<-seq_len(a$max);cdf<-cumsum(k^(-a$exponent))/sum(k^(-a$exponent));vapply(x,function(xx)if(xx<1)0 else if(xx>=a$max)1 else cdf[floor(xx)],numeric(1))},
function(pr,a){k<-seq_len(a$max);cdf<-cumsum(k^(-a$exponent))/sum(k^(-a$exponent));vapply(pr,function(pp)k[which(cdf>=pp)[1L]],numeric(1))},
function(n,a){k<-seq_len(a$max);sample(k,n,replace=TRUE,prob=k^(-a$exponent))}, function(a)c(1,a$max),
function(a){k<-seq_len(a$max);w<-k^(-a$exponent);w<-w/sum(w);m<-sum(k*w);v<-sum((k-m)^2*w);c(mean=m,variance=v,sd=sqrt(v))},
"George Kingsley Zipf popularized an inverse-rank frequency law in linguistics and related empirical systems in the 20th century.",
"Used for rank-frequency and heavy-tailed discrete phenomena. This R4VN implementation is deliberately finite so probabilities, moments, and simulation are well defined for teaching.",
"The finite maximum is part of the model here; changing it changes the normalizing constant and probabilities.", NULL
)
reg$logseries <- S(
"logseries", "Logarithmic series", "discrete", "Discrete - positive counts",
c("log-series", "logarithmic"), list(P("p", "Parameter p", .6, .001, .99, .01, "Shape parameter between 0 and 1.")),
function(x,a) ifelse(x==floor(x)&x>=1,-a$p^x/(x*log(1-a$p)),0),
function(x,a){vapply(x,function(xx){if(xx<1)return(0);k<-seq_len(max(1,floor(xx)));min(1,sum(-a$p^k/(k*log(1-a$p))))},numeric(1))},
function(pr,a){vapply(pr,function(pp){
if (pp <= 0) return(1)
if (pp >= 1) return(Inf)
s<-0;k<-0L
while(s<pp&&k<100000L){k<-k+1L;s<-s-a$p^k/(k*log(1-a$p))}
as.numeric(k)
},numeric(1))},
function(n,a){u<-stats::runif(n);vapply(u,function(pp){s<-0;k<-0L;while(s<pp&&k<100000L){k<-k+1L;s<-s-a$p^k/(k*log(1-a$p))};k},integer(1))},
function(a)c(1,Inf), function(a){m<- -a$p/((1-a$p)*log(1-a$p));v<- -a$p*(a$p+log(1-a$p))/((1-a$p)^2*log(1-a$p)^2);c(mean=m,variance=v,sd=sqrt(v))},
"The logarithmic-series distribution was developed in early-20th-century probability and became well known through R. A. Fisher's work on species abundance.",
"Used for positive, strongly right-skewed counts and historically for species-frequency models.",
"Support starts at 1, unlike Poisson and negative-binomial count distributions which include zero.", NULL
)
reg
}
.r4vn_dist_find <- function(distribution) {
if (missing(distribution) || is.null(distribution) || !length(distribution)) stop("Supply a distribution name, for example `\"normal\"`, `\"binomial\"`, or `\"poisson\"`.", call. = FALSE)
z <- tolower(trimws(as.character(distribution)[1L]))
reg <- .r4vn_dist_registry()
hit <- names(reg)[vapply(reg, function(s) z %in% tolower(s$aliases), logical(1))]
if (!length(hit)) {
labels <- vapply(reg, `[[`, character(1), "id")
stop("Unknown distribution `", distribution, "`. Available names include: ", paste(labels, collapse = ", "), ".", call. = FALSE)
}
reg[[hit[1L]]]
}
.r4vn_dist_validate <- function(spec, a) {
finite1 <- function(nm, positive = FALSE, nonnegative = FALSE, integer = FALSE) {
z <- a[[nm]]
if (!is.numeric(z) || length(z) != 1L || !is.finite(z)) stop("Parameter `", nm, "` must be one finite number.", call. = FALSE)
if (positive && z <= 0) stop("Parameter `", nm, "` must be > 0.", call. = FALSE)
if (nonnegative && z < 0) stop("Parameter `", nm, "` must be >= 0.", call. = FALSE)
if (integer && abs(z-round(z)) > sqrt(.Machine$double.eps)) stop("Parameter `", nm, "` must be an integer.", call. = FALSE)
}
id <- spec$id
positive <- intersect(names(a), c("sd","sdlog","rate","shape","shape1","shape2","scale","df","df1","df2","size","lambda","exponent"))
for (nm in positive) finite1(nm, positive = TRUE)
if ("ncp" %in% names(a)) { finite1("ncp"); if (id %in% c("chisq","f") && a$ncp < 0) stop("`ncp` must be >= 0 for this distribution.", call. = FALSE) }
if ("p" %in% names(a)) { finite1("p"); if (a$p < 0 || a$p > 1 || (id %in% c("geometric","nbinom","logseries") && a$p <= 0) || (id=="logseries" && a$p>=1)) stop("Parameter `p` is outside its valid range.", call. = FALSE) }
if ("zi" %in% names(a)) { finite1("zi"); if (a$zi < 0 || a$zi >= 1) stop("`zi` must satisfy 0 <= zi < 1.", call. = FALSE) }
if (id %in% c("binomial","betabinom")) { finite1("n", nonnegative = TRUE, integer = TRUE); a$n <- as.integer(round(a$n)) }
if (id == "hypergeom") {
for(nm in c("success","failure","draws")) finite1(nm, nonnegative = TRUE, integer = TRUE)
if (a$draws > a$success + a$failure) stop("`draws` cannot exceed the population size.", call. = FALSE)
}
if (id %in% c("uniform","discreteuniform")) {
finite1("min"); finite1("max"); if (a$max <= a$min) stop("`max` must be greater than `min`.", call. = FALSE)
if (id=="discreteuniform" && (a$min!=round(a$min)||a$max!=round(a$max))) stop("Discrete-uniform bounds must be integers.", call. = FALSE)
}
if (id == "triangular") {
finite1("min"); finite1("mode"); finite1("max")
if (!(a$min < a$mode && a$mode < a$max))
stop("Require min < mode < max for the triangular distribution.", call. = FALSE)
}
if (id == "zipf") { finite1("max", positive=TRUE, integer=TRUE); a$max<-as.integer(round(a$max)) }
a
}
.r4vn_dist_args <- function(spec, dots) {
defaults <- setNames(lapply(spec$params, `[[`, "default"), vapply(spec$params, `[[`, character(1), "name"))
if (is.null(names(dots))) names(dots) <- rep("", length(dots))
unnamed <- which(!nzchar(names(dots)))
if (length(unnamed)) stop("Distribution parameters in `...` must be named.", call. = FALSE)
# Friendly aliases are deliberately distribution-specific so that a typo or
# an irrelevant parameter is never silently ignored.
allowed <- names(defaults)
move_alias <- function(old, new) {
if (old %in% names(dots)) {
if (new %in% names(dots))
stop("Supply only one of `", old, "` and `", new, "`.", call. = FALSE)
dots[[new]] <<- dots[[old]]
dots[[old]] <<- NULL
}
}
if ("p" %in% allowed) {
if ("prob" %in% names(dots)) move_alias("prob", "p")
if ("probability" %in% names(dots)) move_alias("probability", "p")
}
if (spec$id %in% c("binomial", "betabinom") && "trials" %in% names(dots))
move_alias("trials", "n")
if (spec$id == "binomial" && "size" %in% names(dots))
move_alias("size", "n")
if (spec$id == "gamma" && "scale" %in% names(dots)) {
if ("rate" %in% names(dots))
stop("For gamma, supply either `rate` or `scale`, not both.", call. = FALSE)
z <- dots$scale
if (!is.numeric(z) || length(z) != 1L || !is.finite(z) || z <= 0)
stop("Gamma `scale` must be one finite number > 0.", call. = FALSE)
dots$rate <- 1 / z
dots$scale <- NULL
}
if (spec$id == "nbinom" && "mu" %in% names(dots)) {
if ("p" %in% names(dots))
stop("For negative binomial, supply either `p` or `mu`, not both.", call. = FALSE)
mu <- dots$mu
size <- dots$size %||% defaults$size
if (!is.numeric(mu) || length(mu) != 1L || !is.finite(mu) || mu < 0)
stop("Negative-binomial `mu` must be one finite number >= 0.", call. = FALSE)
if (!is.numeric(size) || length(size) != 1L || !is.finite(size) || size <= 0)
stop("Negative-binomial `size` must be one finite number > 0.", call. = FALSE)
dots$p <- size / (size + mu)
dots$mu <- NULL
}
if (spec$id == "betabinom" && any(c("p", "rho") %in% names(dots))) {
if (!all(c("p", "rho") %in% names(dots)))
stop("For beta-binomial mean/correlation parameterization, supply both `p` and `rho`.", call. = FALSE)
if (any(c("shape1", "shape2") %in% names(dots)))
stop("For beta-binomial, use either `shape1`/`shape2` or `p`/`rho`, not both.", call. = FALSE)
pp <- dots$p; rho <- dots$rho
if (!is.numeric(pp) || length(pp) != 1L || !is.finite(pp) || pp <= 0 || pp >= 1 ||
!is.numeric(rho) || length(rho) != 1L || !is.finite(rho) || rho <= 0 || rho >= 1)
stop("For beta-binomial `p` and `rho`, require 0 < p < 1 and 0 < rho < 1.", call. = FALSE)
total <- 1 / rho - 1
dots$shape1 <- pp * total
dots$shape2 <- (1 - pp) * total
dots$p <- NULL
dots$rho <- NULL
}
extras <- setdiff(names(dots), allowed)
if (length(extras))
stop("Unused parameter(s) for ", spec$label, ": ", paste(extras, collapse = ", "), ".", call. = FALSE)
for (nm in intersect(names(dots), allowed)) defaults[[nm]] <- dots[[nm]]
.r4vn_dist_validate(spec, defaults)
}
.r4vn_dist_param_text <- function(spec, a, digits = 4L) {
paste(vapply(names(a), function(nm) paste0(nm, " = ", format(signif(a[[nm]], digits), trim=TRUE, scientific=FALSE)), character(1)), collapse = ", ")
}
.r4vn_dist_support_text <- function(spec, a) {
s <- spec$support(a)
lo <- if (is.infinite(s[1])) "-Inf" else format(s[1], trim=TRUE)
hi <- if (is.infinite(s[2])) "Inf" else format(s[2], trim=TRUE)
paste0(if(spec$type=="discrete") "integers in " else "", "[",lo,", ",hi,"]")
}
.r4vn_dist_probability_table <- function(spec, a, x, digits = 6L) {
if (is.null(x) || !length(x)) return(NULL)
x <- as.numeric(x)
rows <- list()
for (xx in x) {
if (spec$type == "discrete") {
eq <- spec$d(xx,a)
lt <- spec$p(ceiling(xx)-1,a)
le <- spec$p(floor(xx),a)
gt <- 1-le
ge <- 1-lt
} else {
eq <- 0
lt <- le <- spec$p(xx,a)
gt <- ge <- 1-le
}
rows[[length(rows)+1L]] <- data.frame(
x = xx,
Quantity = if(spec$type=="discrete") c("P(X = x)","P(X < x)","P(X <= x)","P(X > x)","P(X >= x)") else c("f(x) density","P(X < x)","P(X <= x)","P(X > x)","P(X >= x)"),
Value = c(if(spec$type=="discrete")eq else spec$d(xx,a),lt,le,gt,ge),
stringsAsFactors=FALSE
)
}
out <- do.call(rbind,rows)
# Keep full machine precision in the returned object. `digits` is a display
# control only; rounding here would make programmatic results differ from
# stats::dbinom()/pbinom() and related base-R reference functions.
out$Value <- pmax(0, pmin(ifelse(grepl("density", out$Quantity), Inf, 1), out$Value))
out
}
.r4vn_dist_range_table <- function(spec,a,lower,upper,digits=6L) {
if (is.null(lower) || is.null(upper)) return(NULL)
if (length(lower)!=1L||length(upper)!=1L||!is.finite(lower)||!is.finite(upper)||lower>upper) stop("Supply finite `lower <= upper`.",call.=FALSE)
if(spec$type=="discrete") {
inside <- spec$p(floor(upper),a)-spec$p(ceiling(lower)-1,a)
} else {
inside <- spec$p(upper,a)-spec$p(lower,a)
}
inside <- max(0,min(1,inside))
data.frame(Interval=c("P(lower <= X <= upper)","P(X outside interval)"),
lower=lower,upper=upper,Value=c(inside,1-inside),stringsAsFactors=FALSE)
}
.r4vn_dist_quantile_table <- function(spec,a,probs,digits=6L) {
if (is.null(probs) || !length(probs)) return(NULL)
probs <- as.numeric(probs)
if(any(!is.finite(probs)|probs<0|probs>1))stop("`probs` must contain probabilities between 0 and 1.",call.=FALSE)
data.frame(Probability=probs,Quantile=spec$q(probs,a),stringsAsFactors=FALSE)
}
.r4vn_dist_display_table <- function(x, digits = 6L, columns = NULL) {
if (is.null(x)) return(NULL)
out <- x
if (is.null(columns)) columns <- names(out)[vapply(out, is.numeric, logical(1))]
columns <- intersect(columns, names(out))
for (nm in columns) {
z <- out[[nm]]
ok <- is.finite(z)
z[ok] <- signif(z[ok], digits)
out[[nm]] <- z
}
out
}
.r4vn_dist_grid <- function(spec,a,x=NULL,lower=NULL,upper=NULL,points=501L) {
sup <- spec$support(a)
probs <- if(spec$type=="discrete") c(.001,.999) else c(.005,.995)
qq <- suppressWarnings(spec$q(probs,a))
lo <- if(is.finite(sup[1]))sup[1] else qq[1]
hi <- if(is.finite(sup[2]))sup[2] else qq[2]
extras <- c(x,lower,upper);extras<-extras[is.finite(extras)]
if(length(extras)){lo<-min(lo,extras);hi<-max(hi,extras)}
if(!is.finite(lo)||!is.finite(hi)||lo==hi){lo<- -5;hi<-5}
pad <- .05*(hi-lo); if(spec$type=="continuous"){lo<-if(is.finite(sup[1]))max(sup[1],lo-pad) else lo-pad;hi<-if(is.finite(sup[2]))min(sup[2],hi+pad) else hi+pad}
if(spec$type=="discrete") {
lo <- ceiling(lo); hi <- floor(hi)
if (is.finite(sup[1])) lo <- max(lo, ceiling(sup[1]))
if (is.finite(sup[2])) hi <- min(hi, floor(sup[2]))
if (!is.finite(lo)) lo <- 0
if (!is.finite(hi)) hi <- lo + 50
if (hi < lo) hi <- lo
if(hi-lo>250){center<-suppressWarnings(spec$q(.5,a));if(!is.finite(center))center<-(lo+hi)/2;lo<-max(lo,floor(center-125));hi<-min(hi,ceiling(center+125))}
seq.int(lo,hi)
} else seq(lo,hi,length.out=max(101L,as.integer(points)))
}
#' Explore a probability distribution and calculate probabilities
#'
#' `distdata()` is the command-line probability-distribution calculator for
#' R4VN. It combines distribution properties, point probabilities/densities,
#' tail probabilities, interval probabilities, quantiles, simulation, and a
#' Viewer-ready plot. [distlearn()] provides the interactive Shiny companion.
#'
#' @param distribution Distribution name. Common abbreviations are accepted,
#' for example `"normal"`, `"binomial"`, `"poisson"`, `"t"`, `"chisq"`,
#' `"weibull"`, `"betabinom"`, and `"zip"`.
#' @param ... Named parameters of the selected distribution. For example,
#' `n` and `p` for binomial, `mean` and `sd` for normal, or `lambda` for
#' Poisson. Defaults are teaching-friendly and are shown by `distlearn()`.
#' @param x Optional value(s). For discrete distributions R4VN reports
#' `P(X=x)`, `<`, `<=`, `>`, and `>=`. For continuous distributions it
#' reports density `f(x)` and the corresponding tail probabilities.
#' @param lower,upper Optional interval bounds for an interval probability.
#' @param probs Optional cumulative probabilities for which quantiles are
#' requested, for example `probs = c(.025, .5, .975)`.
#' @param nsim Optional number of random observations to simulate. `0` (the
#' default) does not simulate.
#' @param seed Optional user-supplied random seed for simulation. The default
#' `NULL` does not set a seed.
#' @param plot Logical; retain plot data and display the probability function
#' in the R4VN Viewer.
#' @param digits Number of significant digits in numerical probability output.
#' @param show Logical; open the R4VN Viewer result.
#' @param console Logical; also print the tabular result to the console.
#'
#' @details
#' A particularly useful teaching call is
#' `distdata("binomial", n = 10, x = 3, p = .2)`. With only these inputs R4VN
#' reports `P(X=3)`, `P(X<3)`, `P(X<=3)`, `P(X>3)`, and `P(X>=3)`, together
#' with the distribution's mean, variance, standard deviation, parameters,
#' support, and plot. Supplying `lower` and `upper` adds the probability inside
#' and outside an interval; supplying `probs` adds quantiles.
#'
#' The beta-binomial accepts either `shape1`/`shape2` or the more interpretable
#' pair `p`/`rho`. The negative binomial accepts `size` with either `p` or
#' `mu`. Gamma accepts `rate` or `scale`.
#'
#' @return Invisibly returns an object of classes `r4vn_distdata` and
#' `r4vn_stat`. Its `raw` component contains the distribution specification,
#' parameters, point/range probabilities, quantiles, simulation, and plot
#' grid.
#' @export
#'
#' @examples
#' distdata("binomial", n = 10, x = 3, p = .2, show = FALSE)
#' distdata("binomial", n = 20, p = .35, lower = 5, upper = 10, show = FALSE)
#' distdata("normal", mean = 100, sd = 15, x = 130, show = FALSE)
#' distdata("normal", mean = 100, sd = 15, probs = c(.025, .5, .975), show = FALSE)
#' distdata("poisson", lambda = 2.5, x = 0:4, show = FALSE)
#' distdata("nbinom", size = 2, mu = 5, x = 0:4, show = FALSE)
#' distdata("betabinom", n = 20, p = .3, rho = .1, x = 0:5, show = FALSE)
distdata <- function(distribution, ..., x = NULL, lower = NULL, upper = NULL,
probs = NULL, nsim = 0L, seed = NULL, plot = TRUE,
digits = 6L, show = TRUE, console = FALSE) {
call <- match.call()
if (!is.numeric(digits) || length(digits) != 1L || !is.finite(digits) || digits < 1 || digits > 15)
stop("`digits` must be one integer from 1 to 15.", call. = FALSE)
digits <- as.integer(round(digits))
spec <- .r4vn_dist_find(distribution)
dots <- list(...)
a <- .r4vn_dist_args(spec,dots)
moments <- spec$moments(a)
support <- .r4vn_dist_support_text(spec,a)
params_text <- .r4vn_dist_param_text(spec,a)
prop <- data.frame(
Property = c("Distribution","Type","Parameters","Support","Mean","Variance","Standard deviation"),
Value = c(spec$label,tools::toTitleCase(spec$type),params_text,support,
ifelse(is.na(moments["mean"]),"Not defined / not summarized",format(signif(moments["mean"],digits),trim=TRUE)),
ifelse(is.na(moments["variance"]),"Not defined / not summarized",format(signif(moments["variance"],digits),trim=TRUE)),
ifelse(is.na(moments["sd"]),"Not defined / not summarized",format(signif(moments["sd"],digits),trim=TRUE))),
stringsAsFactors=FALSE
)
point <- .r4vn_dist_probability_table(spec,a,x,digits)
rng <- .r4vn_dist_range_table(spec,a,lower,upper,digits)
quant <- .r4vn_dist_quantile_table(spec,a,probs,digits)
sections <- list("Distribution"=prop)
if(!is.null(point)) sections[["Probability at x"]] <- .r4vn_dist_display_table(point, digits, "Value")
if(!is.null(rng)) sections[["Interval probability"]] <- .r4vn_dist_display_table(rng, digits, "Value")
if(!is.null(quant)) sections[["Quantiles"]] <- .r4vn_dist_display_table(quant, digits, "Quantile")
sim <- NULL
if (!is.numeric(nsim) || length(nsim) != 1L || !is.finite(nsim) || nsim < 0 || abs(nsim - round(nsim)) > sqrt(.Machine$double.eps))
stop("`nsim` must be a non-negative integer.", call. = FALSE)
nsim <- as.integer(round(nsim))
if(nsim>0L){
if(!is.null(seed)){
seed <- .r4vn_seed_normalize(seed)
had<-exists(".Random.seed",envir=.GlobalEnv,inherits=FALSE);if(had)old<-get(".Random.seed",envir=.GlobalEnv,inherits=FALSE)
on.exit({if(had)assign(".Random.seed",old,envir=.GlobalEnv) else if(exists(".Random.seed",envir=.GlobalEnv,inherits=FALSE))rm(".Random.seed",envir=.GlobalEnv)},add=TRUE)
set.seed(seed)
}
sim<-spec$r(nsim,a)
sections[["Simulation"]]<-data.frame(n=nsim,Mean=signif(mean(sim),digits),SD=signif(stats::sd(sim),digits),Minimum=signif(min(sim),digits),Median=signif(stats::median(sim),digits),Maximum=signif(max(sim),digits),stringsAsFactors=FALSE)
}
grid <- if(isTRUE(plot)) .r4vn_dist_grid(spec,a,x,lower,upper) else NULL
raw <- list(spec=spec,parameters=a,moments=moments,point=point,range=rng,quantiles=quant,simulation=sim,grid=grid,x=x,lower=lower,upper=upper)
notes <- c(spec$use, spec$notes)
out <- .r4vn_result(paste0(spec$label," distribution"),sections,notes=notes,raw=raw,call=call)
class(out)<-c("r4vn_distdata",class(out))
.r4vn_show(out,show=show,console=console)
}
#' Plot an R4VN probability-distribution result
#'
#' @param x An object returned by [distdata()].
#' @param type `"density"`/`"mass"` or `"cdf"`.
#' @param main Optional plot title.
#' @param ... Additional graphical parameters passed to base graphics where
#' applicable.
#' @return The input object, invisibly.
#' @export
plot.r4vn_distdata <- function(x, type = c("density","cdf"), main = NULL, ...) {
type <- match.arg(type)
spec<-x$raw$spec;a<-x$raw$parameters;xx<-x$raw$grid
if(is.null(xx))xx<-.r4vn_dist_grid(spec,a,x$raw$x,x$raw$lower,x$raw$upper)
ttl<-main %||% paste0(spec$label,if(type=="cdf")" - cumulative distribution" else if(spec$type=="discrete")" - probability mass" else " - density")
if(type=="cdf"){
yy<-spec$p(xx,a)
if(spec$type=="discrete")graphics::plot(xx,yy,type="s",xlab="x",ylab="F(x)",main=ttl,...) else graphics::plot(xx,yy,type="l",xlab="x",ylab="F(x)",main=ttl,...)
graphics::abline(h=c(0,1),lty=3)
} else if(spec$type=="discrete"){
yy<-spec$d(xx,a);graphics::plot(xx,yy,type="h",lwd=3,xlab="x",ylab="P(X = x)",main=ttl,...);graphics::points(xx,yy,pch=16)
} else {
yy<-spec$d(xx,a);graphics::plot(xx,yy,type="l",lwd=2,xlab="x",ylab="Density",main=ttl,...)
}
if(!is.null(x$raw$x)){
xv<-x$raw$x[is.finite(x$raw$x)];if(length(xv))graphics::abline(v=xv,lty=2)
}
if(!is.null(x$raw$lower)&&!is.null(x$raw$upper))graphics::abline(v=c(x$raw$lower,x$raw$upper),lty=2)
invisible(x)
}
.r4vn_dist_svg <- function(x,type="density") {
f<-tempfile(fileext=".svg")
grDevices::svg(f,width=7.5,height=4.3,pointsize=11)
on.exit({try(grDevices::dev.off(),silent=TRUE);unlink(f)},add=TRUE)
plot(x,type=type)
grDevices::dev.off()
txt<-readLines(f,warn=FALSE,encoding="UTF-8")
start<-grep("<svg",txt,fixed=TRUE)[1L]
if(is.na(start))return("")
paste(txt[start:length(txt)],collapse="\n")
}
.r4vn_viewer_distdata <- function(x) {
spec<-x$raw$spec
hero<-paste0('<section class="r4vn-section"><div style="display:grid;grid-template-columns:1fr 1fr;gap:12px">',
'<div><h2>Why this distribution?</h2><p>',.r4vn_view_escape(spec$use),'</p></div>',
'<div><h2>Historical note</h2><p>',.r4vn_view_escape(spec$history),'</p></div></div></section>')
graph<-if(is.null(x$raw$grid))"" else paste0('<section class="r4vn-section"><h2>Distribution plot</h2>',.r4vn_dist_svg(x,"density"),'</section>')
blocks<-paste0(vapply(names(x$sections),function(nm).r4vn_view_section(nm,x$sections[[nm]]),character(1)),collapse="")
formula<-if(is.null(spec$formula)||!nzchar(spec$formula))"" else paste0('<section class="r4vn-section"><h2>Formula / definition</h2><pre>',.r4vn_view_escape(spec$formula),'</pre></section>')
.r4vn_view_document(x$title,paste0(hero,graph,blocks,formula),notes=x$notes,subtitle="R4VN distribution calculator",prefix="r4vn-distdata-")
}
Try the R4VN package in your browser
Any scripts or data that you put into this service are public.
R4VN documentation built on Sept. 30, 2026, 5:13 p.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.
Defines functions .r4vn_viewer_distdata .r4vn_dist_svg plot.r4vn_distdata distdata .r4vn_dist_grid .r4vn_dist_display_table .r4vn_dist_quantile_table .r4vn_dist_range_table .r4vn_dist_probability_table .r4vn_dist_support_text .r4vn_dist_param_text .r4vn_dist_args .r4vn_dist_validate .r4vn_dist_find .r4vn_dist_registry .r4vn_dist_safe_var .r4vn_dist_spec .r4vn_dist_param
Documented in distdata plot.r4vn_distdata
# ============================================================================
# R4VN probability distributions: command engine and teaching calculator
# ============================================================================
.r4vn_dist_param <- function(name, label, default, min = NULL, max = NULL, step = NULL, description = "") {
list(name = name, label = label, default = default, min = min, max = max, step = step, description = description)
}
.r4vn_dist_spec <- function(id, label, type, category, aliases, params,
d, p, q, r, support, moments, history, use, notes = NULL,
formula = NULL) {
list(id = id, label = label, type = type, category = category,
aliases = unique(c(id, aliases)), params = params,
d = d, p = p, q = q, r = r, support = support, moments = moments,
history = history, use = use, notes = notes, formula = formula)
}
.r4vn_dist_safe_var <- function(x) if (is.finite(x) && x >= 0) x else NA_real_
.r4vn_dist_registry <- function() {
P <- .r4vn_dist_param
S <- .r4vn_dist_spec
reg <- list()
reg$normal <- S(
"normal", "Normal (Gaussian)", "continuous", "Continuous - symmetric",
c("norm", "gaussian"),
list(P("mean", "Mean (mu)", 0, -10, 10, .1, "Center of the distribution."),
P("sd", "Standard deviation (sigma)", 1, .05, 10, .05, "Controls spread; must be positive.")),
function(x,a) stats::dnorm(x,a$mean,a$sd), function(x,a) stats::pnorm(x,a$mean,a$sd),
function(pr,a) stats::qnorm(pr,a$mean,a$sd), function(n,a) stats::rnorm(n,a$mean,a$sd),
function(a) c(-Inf,Inf), function(a) c(mean=a$mean,variance=a$sd^2,sd=a$sd),
"The bell-shaped law emerged from work by Abraham de Moivre on the binomial approximation in the 18th century and was developed further by Laplace and Gauss in probability and error theory.",
"Common for approximately symmetric continuous measurements, sampling distributions, measurement error, and many regression assumptions. It should not be chosen only because a variable is continuous.",
"About 68%, 95%, and 99.7% of values lie within 1, 2, and 3 standard deviations of the mean when the distribution is exactly normal.",
"f(x) = exp[-(x-mu)^2/(2 sigma^2)] / (sigma sqrt(2 pi))"
)
reg$lognormal <- S(
"lognormal", "Log-normal", "continuous", "Continuous - positive",
c("lnorm", "log-normal"),
list(P("meanlog", "Mean on log scale", 0, -5, 5, .1, "Mean of log(X)."),
P("sdlog", "SD on log scale", 1, .05, 3, .05, "Standard deviation of log(X).")),
function(x,a) stats::dlnorm(x,a$meanlog,a$sdlog), function(x,a) stats::plnorm(x,a$meanlog,a$sdlog),
function(pr,a) stats::qlnorm(pr,a$meanlog,a$sdlog), function(n,a) stats::rlnorm(n,a$meanlog,a$sdlog),
function(a)c(0,Inf), function(a){m<-exp(a$meanlog+a$sdlog^2/2);v<-(exp(a$sdlog^2)-1)*exp(2*a$meanlog+a$sdlog^2);c(mean=m,variance=v,sd=sqrt(v))},
"The log-normal family became important through 19th- and early-20th-century work on multiplicative variation; it describes a variable whose logarithm is normally distributed.",
"Useful for positive right-skewed quantities produced by multiplicative processes, such as some biological measurements, costs, incubation times, and environmental concentrations.",
"The parameters describe log(X), not X. The median on the original scale is exp(meanlog).",
"log(X) ~ Normal(meanlog, sdlog^2)"
)
reg$uniform <- S(
"uniform", "Uniform", "continuous", "Continuous - bounded",
c("unif"),
list(P("min", "Minimum (a)", 0, -10, 10, .1, "Lower support boundary."), P("max", "Maximum (b)", 1, -10, 20, .1, "Upper support boundary; must exceed a.")),
function(x,a) stats::dunif(x,a$min,a$max), function(x,a) stats::punif(x,a$min,a$max), function(pr,a) stats::qunif(pr,a$min,a$max), function(n,a) stats::runif(n,a$min,a$max),
function(a)c(a$min,a$max), function(a)c(mean=(a$min+a$max)/2,variance=(a$max-a$min)^2/12,sd=(a$max-a$min)/sqrt(12)),
"Uniform models are among the oldest probability models and formalize equal density over a bounded interval.",
"Useful as a simple bounded model, for random-number generation, simulation, and as a non-informative building block. It is usually unrealistic for biological measurements unless equal density is scientifically plausible.",
"Every interval of the same width inside [a,b] has the same probability.", "f(x)=1/(b-a), a <= x <= b"
)
reg$exponential <- S(
"exponential", "Exponential", "continuous", "Continuous - positive",
c("exp"), list(P("rate", "Rate (lambda)", 1, .01, 10, .05, "Event rate per unit time; mean equals 1/rate.")),
function(x,a) stats::dexp(x,a$rate), function(x,a) stats::pexp(x,a$rate), function(pr,a) stats::qexp(pr,a$rate), function(n,a) stats::rexp(n,a$rate),
function(a)c(0,Inf), function(a)c(mean=1/a$rate,variance=1/a$rate^2,sd=1/a$rate),
"The exponential distribution arose naturally from Poisson-process theory and waiting-time problems.",
"Used for positive waiting times when the instantaneous hazard is constant. In survival analysis, that constant-hazard assumption should be considered explicitly.",
"It is memoryless: conditional remaining time has the same distribution regardless of elapsed time.", "f(x)=lambda exp(-lambda x), x >= 0"
)
reg$gamma <- S(
"gamma", "Gamma", "continuous", "Continuous - positive",
c("gam"), list(P("shape", "Shape (alpha)", 2, .05, 15, .05, "Controls skewness and shape."), P("rate", "Rate (beta)", 1, .01, 10, .05, "Inverse scale; mean = shape/rate.")),
function(x,a) stats::dgamma(x,a$shape,rate=a$rate), function(x,a) stats::pgamma(x,a$shape,rate=a$rate), function(pr,a) stats::qgamma(pr,a$shape,rate=a$rate), function(n,a) stats::rgamma(n,a$shape,rate=a$rate),
function(a)c(0,Inf), function(a)c(mean=a$shape/a$rate,variance=a$shape/a$rate^2,sd=sqrt(a$shape)/a$rate),
"The family is tied to Euler's gamma function and became central in waiting-time and Pearson-family models.",
"Useful for positive right-skewed continuous outcomes, waiting times to multiple events, costs, and as a flexible GLM distribution with a log link.",
"R4VN uses shape-rate parameterization. `scale` can be supplied to distdata() as an alternative, where rate = 1/scale.", "f(x)=beta^alpha x^(alpha-1) exp(-beta x) / Gamma(alpha)"
)
reg$beta <- S(
"beta", "Beta", "continuous", "Continuous - bounded",
c("beta"), list(P("shape1", "Shape 1 (alpha)", 2, .05, 15, .05, "Controls behavior near 0 and overall shape."), P("shape2", "Shape 2 (beta)", 2, .05, 15, .05, "Controls behavior near 1 and overall shape.")),
function(x,a) stats::dbeta(x,a$shape1,a$shape2), function(x,a) stats::pbeta(x,a$shape1,a$shape2), function(pr,a) stats::qbeta(pr,a$shape1,a$shape2), function(n,a) stats::rbeta(n,a$shape1,a$shape2),
function(a)c(0,1), function(a){m<-a$shape1/(a$shape1+a$shape2);v<-a$shape1*a$shape2/((a$shape1+a$shape2)^2*(a$shape1+a$shape2+1));c(mean=m,variance=v,sd=sqrt(v))},
"The beta family is based on Euler's beta function and later became a core bounded distribution in probability and Bayesian statistics.",
"Useful for proportions or probabilities that vary continuously between 0 and 1 and as a conjugate prior for a Bernoulli/binomial probability.",
"alpha=beta gives symmetry; values below 1 can create U-shaped densities, while values above 1 often create unimodal densities.", "f(x)=x^(alpha-1)(1-x)^(beta-1)/B(alpha,beta)"
)
reg$chisq <- S(
"chisq", "Chi-square", "continuous", "Continuous - sampling distributions",
c("chi-square", "chi square", "chi2"), list(P("df", "Degrees of freedom", 4, .1, 50, .1, "Controls shape and spread."), P("ncp", "Noncentrality", 0, 0, 25, .1, "0 gives the central chi-square distribution.")),
function(x,a) stats::dchisq(x,a$df,ncp=a$ncp), function(x,a) stats::pchisq(x,a$df,ncp=a$ncp), function(pr,a) stats::qchisq(pr,a$df,ncp=a$ncp), function(n,a) stats::rchisq(n,a$df,ncp=a$ncp),
function(a)c(0,Inf), function(a){if(a$ncp==0)c(mean=a$df,variance=2*a$df,sd=sqrt(2*a$df)) else c(mean=a$df+a$ncp,variance=2*(a$df+2*a$ncp),sd=sqrt(2*(a$df+2*a$ncp)))},
"Karl Pearson introduced the chi-square goodness-of-fit framework around 1900; the distribution also arises as a sum of squared standard normal variables.",
"Used in categorical-data tests, variance inference under normality, likelihood theory, and many asymptotic test statistics.",
"For the central distribution, the mean is df and variance is 2*df. A noncentrality parameter is useful in power calculations.", "X = sum(Z_i^2) for df independent standard normal variables"
)
reg$t <- S(
"t", "Student's t", "continuous", "Continuous - sampling distributions",
c("student", "student-t", "student t"), list(P("df", "Degrees of freedom", 5, .5, 100, .5, "Smaller df gives heavier tails."), P("ncp", "Noncentrality", 0, -10, 10, .1, "0 gives the central t distribution.")),
function(x,a) stats::dt(x,a$df,ncp=a$ncp), function(x,a) stats::pt(x,a$df,ncp=a$ncp), function(pr,a) stats::qt(pr,a$df,ncp=a$ncp), function(n,a) stats::rt(n,a$df,ncp=a$ncp),
function(a)c(-Inf,Inf), function(a){if(a$ncp!=0)return(c(mean=NA_real_,variance=NA_real_,sd=NA_real_));m<-if(a$df>1)0 else NA_real_;v<-if(a$df>2)a$df/(a$df-2) else if(a$df>1) Inf else NA_real_;c(mean=m,variance=v,sd=sqrt(v))},
"William Sealy Gosset published the distribution under the pseudonym 'Student' in 1908 while working at Guinness; Fisher later developed its role in classical inference.",
"Used for inference on means and regression coefficients when a standard error is estimated, especially in small samples under normal-model assumptions.",
"As df increases, the central t distribution approaches the standard normal; small df produce heavier tails.", NULL
)
reg$f <- S(
"f", "F", "continuous", "Continuous - sampling distributions",
c("f-distribution", "f distribution"), list(P("df1", "Numerator df", 5, .5, 50, .5, "Numerator degrees of freedom."), P("df2", "Denominator df", 20, .5, 100, .5, "Denominator degrees of freedom."), P("ncp", "Noncentrality", 0, 0, 25, .1, "0 gives the central F distribution.")),
function(x,a) stats::df(x,a$df1,a$df2,ncp=a$ncp), function(x,a) stats::pf(x,a$df1,a$df2,ncp=a$ncp), function(pr,a) stats::qf(pr,a$df1,a$df2,ncp=a$ncp), function(n,a) stats::rf(n,a$df1,a$df2,ncp=a$ncp),
function(a)c(0,Inf), function(a){if(a$ncp!=0)return(c(mean=NA_real_,variance=NA_real_,sd=NA_real_));m<-if(a$df2>2)a$df2/(a$df2-2) else NA_real_;v<-if(a$df2>4)2*a$df2^2*(a$df1+a$df2-2)/(a$df1*(a$df2-2)^2*(a$df2-4)) else NA_real_;c(mean=m,variance=v,sd=sqrt(v))},
"The F distribution is associated with Ronald A. Fisher's variance-ratio methods and was named F by George W. Snedecor.",
"Used in ANOVA, nested-model tests, variance-ratio problems, and regression/model comparison. Noncentral F distributions are common in power calculations.",
"It is a ratio of two appropriately scaled independent chi-square variables in the central case.", NULL
)
reg$weibull <- S(
"weibull", "Weibull", "continuous", "Continuous - positive",
c("weib"), list(P("shape", "Shape", 1.5, .1, 8, .05, "Determines how hazard changes over time."), P("scale", "Scale", 1, .05, 10, .05, "Sets the time/measurement scale.")),
function(x,a) stats::dweibull(x,a$shape,a$scale), function(x,a) stats::pweibull(x,a$shape,a$scale), function(pr,a) stats::qweibull(pr,a$shape,a$scale), function(n,a) stats::rweibull(n,a$shape,a$scale),
function(a)c(0,Inf), function(a){m<-a$scale*gamma(1+1/a$shape);v<-a$scale^2*(gamma(1+2/a$shape)-gamma(1+1/a$shape)^2);c(mean=m,variance=v,sd=sqrt(v))},
"The distribution is named for Waloddi Weibull, who demonstrated its broad applicability in reliability work in the 20th century.",
"Widely used for survival and reliability data. Shape <1 gives decreasing hazard, shape=1 gives the exponential model, and shape >1 gives increasing hazard.",
"Its flexible hazard pattern makes it useful when constant hazard is implausible but monotonic hazard is reasonable.", NULL
)
reg$logistic <- S(
"logistic", "Logistic", "continuous", "Continuous - symmetric",
c("logis"), list(P("location", "Location", 0, -10, 10, .1, "Center and median."), P("scale", "Scale", 1, .05, 10, .05, "Positive spread parameter.")),
function(x,a) stats::dlogis(x,a$location,a$scale), function(x,a) stats::plogis(x,a$location,a$scale), function(pr,a) stats::qlogis(pr,a$location,a$scale), function(n,a) stats::rlogis(n,a$location,a$scale),
function(a)c(-Inf,Inf), function(a)c(mean=a$location,variance=pi^2*a$scale^2/3,sd=pi*a$scale/sqrt(3)),
"The logistic curve is associated with Pierre-Francois Verhulst's 19th-century population-growth work; the corresponding probability distribution later became important in statistical modeling.",
"Useful as a symmetric heavy-tailed alternative to the normal and fundamental to logistic regression through its cumulative distribution/logit relationship.",
"Compared with a normal distribution of similar spread, the logistic has heavier tails.", NULL
)
reg$cauchy <- S(
"cauchy", "Cauchy", "continuous", "Continuous - symmetric",
c("lorentz"), list(P("location", "Location", 0, -10, 10, .1, "Median and center of symmetry."), P("scale", "Scale", 1, .05, 10, .05, "Positive scale parameter.")),
function(x,a) stats::dcauchy(x,a$location,a$scale), function(x,a) stats::pcauchy(x,a$location,a$scale), function(pr,a) stats::qcauchy(pr,a$location,a$scale), function(n,a) stats::rcauchy(n,a$location,a$scale),
function(a)c(-Inf,Inf), function(a)c(mean=NA_real_,variance=NA_real_,sd=NA_real_),
"Named for Augustin-Louis Cauchy, this heavy-tailed law also appears as the ratio of two independent standard normal variables.",
"Useful as a heavy-tailed model and in robust/Bayesian work, but not as a routine substitute for normality. Its extreme tails materially affect summaries and inference.",
"The ordinary mean and variance do not exist, even though the distribution is symmetric around its location parameter.", NULL
)
reg$laplace <- S(
"laplace", "Laplace (double exponential)", "continuous", "Continuous - symmetric",
c("double exponential"), list(P("location", "Location", 0, -10, 10, .1, "Center and median."), P("scale", "Scale", 1, .05, 10, .05, "Mean absolute deviation scale; must be positive.")),
function(x,a) exp(-abs(x-a$location)/a$scale)/(2*a$scale),
function(x,a) ifelse(x<a$location,.5*exp((x-a$location)/a$scale),1-.5*exp(-(x-a$location)/a$scale)),
function(pr,a) a$location + a$scale*ifelse(pr<.5,log(2*pr),-log(2*(1-pr))),
function(n,a){u<-stats::runif(n);a$location+a$scale*ifelse(u<.5,log(2*u),-log(2*(1-u)))},
function(a)c(-Inf,Inf), function(a)c(mean=a$location,variance=2*a$scale^2,sd=sqrt(2)*a$scale),
"The distribution is associated with Pierre-Simon Laplace and is also called the double-exponential distribution.",
"Useful for symmetric data with a sharp center and heavier tails than the normal; it also underlies least-absolute-deviation ideas and the lasso's common Bayesian prior interpretation.",
"The scale parameter is not the standard deviation; SD = sqrt(2)*scale.", NULL
)
reg$gumbel <- S(
"gumbel", "Gumbel (extreme value type I)", "continuous", "Continuous - extreme values",
c("ev1", "extreme value"), list(P("location", "Location", 0, -10, 10, .1, "Shifts the distribution."), P("scale", "Scale", 1, .05, 10, .05, "Positive spread parameter.")),
function(x,a){z<-(x-a$location)/a$scale;exp(-(z+exp(-z)))/a$scale},
function(x,a){z<-(x-a$location)/a$scale;exp(-exp(-z))},
function(pr,a)a$location-a$scale*log(-log(pr)), function(n,a){u<-stats::runif(n);a$location-a$scale*log(-log(u))},
function(a)c(-Inf,Inf), function(a){eg<-0.5772156649015329;m<-a$location+eg*a$scale;v<-pi^2*a$scale^2/6;c(mean=m,variance=v,sd=sqrt(v))},
"Named for Emil Gumbel, a major developer of extreme-value statistics in the 20th century.",
"Used to model maxima or minima (after an appropriate sign convention), such as extreme environmental measurements or maximum loads.",
"This implementation uses the distribution for maxima with CDF exp(-exp(-(x-location)/scale)).", NULL
)
reg$pareto <- S(
"pareto", "Pareto type I", "continuous", "Continuous - heavy tailed",
c("pareto1"), list(P("shape", "Shape (alpha)", 3, .1, 10, .05, "Tail index; larger values give lighter tails."), P("scale", "Minimum (x_m)", 1, .05, 10, .05, "Lower support boundary.")),
function(x,a){out<-rep(0,length(x));ok<-is.finite(x)&x>=a$scale;out[ok]<-a$shape*a$scale^a$shape/x[ok]^(a$shape+1);out[is.infinite(x)&x>0]<-0;out},
function(x,a){out<-rep(0,length(x));ok<-is.finite(x)&x>=a$scale;out[ok]<-1-(a$scale/x[ok])^a$shape;out[is.infinite(x)&x>0]<-1;out},
function(pr,a)a$scale/(1-pr)^(1/a$shape), function(n,a){u<-stats::runif(n);a$scale/(1-u)^(1/a$shape)},
function(a)c(a$scale,Inf), function(a){m<-if(a$shape>1)a$shape*a$scale/(a$shape-1) else Inf;v<-if(a$shape>2)a$shape*a$scale^2/((a$shape-1)^2*(a$shape-2)) else Inf;c(mean=m,variance=v,sd=sqrt(v))},
"Vilfredo Pareto used power-law ideas in studying wealth and income; the Pareto family later became a standard heavy-tail model.",
"Useful for strongly right-skewed heavy-tailed quantities and threshold exceedances in some contexts. Tail assumptions should be checked because moments may not exist.",
"The mean is infinite when shape <=1 and the variance is infinite when shape <=2.", NULL
)
reg$rayleigh <- S(
"rayleigh", "Rayleigh", "continuous", "Continuous - positive",
c("rayleigh"), list(P("scale", "Scale (sigma)", 1, .05, 10, .05, "Positive scale parameter.")),
function(x,a) ifelse(x<0,0,x/a$scale^2*exp(-x^2/(2*a$scale^2))),
function(x,a) ifelse(x<0,0,1-exp(-x^2/(2*a$scale^2))),
function(pr,a)a$scale*sqrt(-2*log(1-pr)), function(n,a)a$scale*sqrt(-2*log(stats::runif(n))),
function(a)c(0,Inf), function(a){m<-a$scale*sqrt(pi/2);v<-(2-pi/2)*a$scale^2;c(mean=m,variance=v,sd=sqrt(v))},
"Named for Lord Rayleigh; the distribution arises as the magnitude of a two-dimensional vector with independent zero-mean normal components of equal variance.",
"Used for magnitudes, signal amplitudes, wind-speed-like models, and radial error when two orthogonal components are Gaussian.",
"It is positive and right-skewed; the single scale parameter controls both center and spread.", NULL
)
reg$triangular <- S(
"triangular", "Triangular", "continuous", "Continuous - bounded",
c("triangle"), list(P("min", "Minimum (a)", 0, -10, 10, .1, "Lower bound."), P("mode", "Mode (c)", .5, -10, 20, .1, "Most likely value; must lie between min and max."), P("max", "Maximum (b)", 1, -10, 20, .1, "Upper bound.")),
function(x,a){out<-rep(0,length(x));i<-x>=a$min&x<=a$mode;out[i]<-2*(x[i]-a$min)/((a$max-a$min)*(a$mode-a$min));j<-x>a$mode&x<=a$max;out[j]<-2*(a$max-x[j])/((a$max-a$min)*(a$max-a$mode));out},
function(x,a){out<-ifelse(x<a$min,0,ifelse(x>=a$max,1,NA_real_));i<-x>=a$min&x<=a$mode;out[i]<-(x[i]-a$min)^2/((a$max-a$min)*(a$mode-a$min));j<-x>a$mode&x<a$max;out[j]<-1-(a$max-x[j])^2/((a$max-a$min)*(a$max-a$mode));out},
function(pr,a){fc<-(a$mode-a$min)/(a$max-a$min);ifelse(pr<fc,a$min+sqrt(pr*(a$max-a$min)*(a$mode-a$min)),a$max-sqrt((1-pr)*(a$max-a$min)*(a$max-a$mode)))},
function(n,a){u<-stats::runif(n);fc<-(a$mode-a$min)/(a$max-a$min);ifelse(u<fc,a$min+sqrt(u*(a$max-a$min)*(a$mode-a$min)),a$max-sqrt((1-u)*(a$max-a$min)*(a$max-a$mode)))},
function(a)c(a$min,a$max), function(a){m<-(a$min+a$mode+a$max)/3;v<-(a$min^2+a$mode^2+a$max^2-a$min*a$mode-a$min*a$max-a$mode*a$max)/18;c(mean=m,variance=v,sd=sqrt(v))},
"The triangular distribution is a simple bounded model long used in simulation and project-risk settings when only a minimum, mode, and maximum are elicited.",
"Useful for transparent teaching and rough uncertainty modeling when information is limited to plausible bounds and a most likely value; not a substitute for empirical distribution fitting.",
"The mode must lie between the lower and upper bounds.", NULL
)
reg$kumaraswamy <- S(
"kumaraswamy", "Kumaraswamy", "continuous", "Continuous - bounded",
c("kuma"), list(P("shape1", "Shape a", 2, .05, 15, .05, "First shape parameter."), P("shape2", "Shape b", 2, .05, 15, .05, "Second shape parameter.")),
function(x,a) ifelse(x<0|x>1,0,a$shape1*a$shape2*x^(a$shape1-1)*(1-x^a$shape1)^(a$shape2-1)),
function(x,a) ifelse(x<=0,0,ifelse(x>=1,1,1-(1-x^a$shape1)^a$shape2)),
function(pr,a)(1-(1-pr)^(1/a$shape2))^(1/a$shape1), function(n,a){u<-stats::runif(n);(1-(1-u)^(1/a$shape2))^(1/a$shape1)},
function(a)c(0,1), function(a){m<-a$shape2*beta(1+1/a$shape1,a$shape2);e2<-a$shape2*beta(1+2/a$shape1,a$shape2);v<-e2-m^2;c(mean=m,variance=v,sd=sqrt(max(0,v)))},
"P. Kumaraswamy proposed this bounded two-shape distribution in 1980 for hydrologic applications.",
"Useful as a flexible beta-like model on (0,1), particularly in simulation because its CDF and quantile function have simple closed forms.",
"It resembles the beta family in flexibility but uses different shape parameters and is not interchangeable with beta regression parameterizations.", NULL
)
reg$loglogistic <- S(
"loglogistic", "Log-logistic", "continuous", "Continuous - positive",
c("log-logistic", "fisk"), list(P("shape", "Shape", 2.5, .1, 10, .05, "Tail/hazard shape."), P("scale", "Scale", 1, .05, 10, .05, "Median and positive scale.")),
function(x,a){out<-rep(0,length(x));ok<-is.finite(x)&x>0;z<-x[ok]/a$scale;out[ok]<-(a$shape/a$scale)*z^(a$shape-1)/(1+z^a$shape)^2;out[is.infinite(x)&x>0]<-0;out},
function(x,a){out<-rep(0,length(x));ok<-is.finite(x)&x>0;z<-x[ok]/a$scale;out[ok]<-1/(1+z^(-a$shape));out[is.infinite(x)&x>0]<-1;out},
function(pr,a)a$scale*(pr/(1-pr))^(1/a$shape), function(n,a){u<-stats::runif(n);a$scale*(u/(1-u))^(1/a$shape)},
function(a)c(0,Inf), function(a){m<-if(a$shape>1)a$scale*(pi/a$shape)/sin(pi/a$shape) else Inf;e2<-if(a$shape>2)a$scale^2*(2*pi/a$shape)/sin(2*pi/a$shape) else Inf;v<-if(is.finite(e2))e2-m^2 else Inf;c(mean=m,variance=v,sd=sqrt(v))},
"The log-logistic distribution results when log(X) follows a logistic distribution and is also known as the Fisk distribution in some fields.",
"Used for positive duration/survival data when a non-monotone hazard may be plausible, and for some income or size distributions.",
"The scale parameter is the median. The mean exists only for shape >1 and variance only for shape >2.", NULL
)
reg$halfnormal <- S(
"halfnormal", "Half-normal", "continuous", "Continuous - positive",
c("half-normal"), list(P("sd", "Underlying SD (sigma)", 1, .05, 10, .05, "SD of the zero-mean normal variable before taking absolute value.")),
function(x,a) ifelse(x<0,0,2*stats::dnorm(x,0,a$sd)), function(x,a) ifelse(x<0,0,2*stats::pnorm(x/a$sd)-1),
function(pr,a)a$sd*stats::qnorm((pr+1)/2), function(n,a)abs(stats::rnorm(n,0,a$sd)), function(a)c(0,Inf),
function(a){m<-a$sd*sqrt(2/pi);v<-a$sd^2*(1-2/pi);c(mean=m,variance=v,sd=sqrt(v))},
"The half-normal is the distribution of the absolute value of a zero-mean normal variable.",
"Useful for non-negative magnitudes, random-effect scale priors, measurement-error magnitudes, and teaching transformations of normal variables.",
"Its parameter is the SD of the underlying normal, not the SD of the half-normal itself.", NULL
)
# Discrete distributions ---------------------------------------------------
reg$bernoulli <- S(
"bernoulli", "Bernoulli", "discrete", "Discrete - counts/events",
c("bern"), list(P("p", "Probability of success", .5, 0, 1, .01, "Probability that X=1.")),
function(x,a) stats::dbinom(x,1,a$p), function(x,a) stats::pbinom(x,1,a$p), function(pr,a) stats::qbinom(pr,1,a$p), function(n,a) stats::rbinom(n,1,a$p),
function(a)c(0,1), function(a)c(mean=a$p,variance=a$p*(1-a$p),sd=sqrt(a$p*(1-a$p))),
"Named for Jacob Bernoulli, whose work on repeated binary trials was published posthumously in Ars Conjectandi in 1713.",
"Models one binary trial: event/no event, success/failure, diseased/not diseased. Repeated independent Bernoulli trials lead to the binomial distribution.",
"X takes only 0 and 1; p is both P(X=1) and E(X).", "P(X=x)=p^x(1-p)^(1-x), x in {0,1}"
)
reg$binomial <- S(
"binomial", "Binomial", "discrete", "Discrete - counts/events",
c("binom"), list(P("n", "Number of trials (n)", 10, 1, 200, 1, "Number of independent Bernoulli trials."), P("p", "Probability of success", .5, 0, 1, .01, "Success probability for each trial.")),
function(x,a) stats::dbinom(x,a$n,a$p), function(x,a) stats::pbinom(x,a$n,a$p), function(pr,a) stats::qbinom(pr,a$n,a$p), function(nn,a) stats::rbinom(nn,a$n,a$p),
function(a)c(0,a$n), function(a)c(mean=a$n*a$p,variance=a$n*a$p*(1-a$p),sd=sqrt(a$n*a$p*(1-a$p))),
"The binomial law grew from early probability work on repeated trials, especially Jacob Bernoulli's law of large numbers and binomial-trial reasoning.",
"Used for the number of events in a fixed number of independent trials when each trial has the same event probability. Examples include positive tests among n specimens or responders among n independent participants under a simple model.",
"Key assumptions are fixed n, two outcomes per trial, constant p, and independence. Overdispersion or clustering can invalidate the simple binomial model.", "P(X=x)=choose(n,x) p^x (1-p)^(n-x)"
)
reg$poisson <- S(
"poisson", "Poisson", "discrete", "Discrete - counts/events",
c("pois"), list(P("lambda", "Mean/rate (lambda)", 3, .01, 30, .05, "Expected count in the specified exposure interval.")),
function(x,a) stats::dpois(x,a$lambda), function(x,a) stats::ppois(x,a$lambda), function(pr,a) stats::qpois(pr,a$lambda), function(n,a) stats::rpois(n,a$lambda),
function(a)c(0,Inf), function(a)c(mean=a$lambda,variance=a$lambda,sd=sqrt(a$lambda)),
"Simeon Denis Poisson described the distribution in 1837 as a limiting law for rare events.",
"Used for counts over time, person-time, area, or another exposure when events occur independently at a stable rate. Poisson regression extends this idea to rate modeling with covariates and offsets.",
"The basic Poisson distribution has mean equal to variance; extra-Poisson variation suggests overdispersion, clustering, heterogeneity, or another count model.", "P(X=x)=exp(-lambda) lambda^x/x!"
)
reg$geometric <- S(
"geometric", "Geometric", "discrete", "Discrete - waiting counts",
c("geom"), list(P("p", "Success probability", .25, .001, 1, .01, "Probability of success on each independent trial.")),
function(x,a) stats::dgeom(x,a$p), function(x,a) stats::pgeom(x,a$p), function(pr,a) stats::qgeom(pr,a$p), function(n,a) stats::rgeom(n,a$p),
function(a)c(0,Inf), function(a)c(mean=(1-a$p)/a$p,variance=(1-a$p)/a$p^2,sd=sqrt(1-a$p)/a$p),
"The geometric distribution is one of the classical waiting-time distributions arising from Bernoulli trials.",
"Used for the number of failures before the first success under independent trials with constant success probability.",
"R4VN follows R's convention: X counts failures before the first success, so support starts at 0. Some textbooks count the trial number instead and therefore start at 1.", NULL
)
reg$nbinom <- S(
"nbinom", "Negative binomial", "discrete", "Discrete - counts/events",
c("negative binomial", "negbin"), list(P("size", "Size/dispersion", 2, .05, 30, .05, "Positive shape/dispersion parameter."), P("p", "Success probability", .4, .001, .999, .01, "R probability parameter; mean=size*(1-p)/p.")),
function(x,a) stats::dnbinom(x,size=a$size,prob=a$p), function(x,a) stats::pnbinom(x,size=a$size,prob=a$p), function(pr,a) stats::qnbinom(pr,size=a$size,prob=a$p), function(n,a) stats::rnbinom(n,size=a$size,prob=a$p),
function(a)c(0,Inf), function(a){m<-a$size*(1-a$p)/a$p;v<-a$size*(1-a$p)/a$p^2;c(mean=m,variance=v,sd=sqrt(v))},
"The negative-binomial family developed from classical repeated-trial problems and later became important as a Poisson-gamma mixture for overdispersed counts.",
"Used for overdispersed count data and heterogeneous event rates. Negative-binomial regression is a common alternative when Poisson variance is too small.",
"distdata() also accepts `mu` instead of `p`; it converts using p=size/(size+mu).", NULL
)
reg$hypergeom <- S(
"hypergeom", "Hypergeometric", "discrete", "Discrete - sampling without replacement",
c("hypergeometric", "hyper"), list(P("success", "Successes in population", 10, 0, 100, 1, "Number of population items classified as success."), P("failure", "Failures in population", 20, 0, 100, 1, "Number of population items classified as failure."), P("draws", "Sample size", 5, 0, 100, 1, "Number drawn without replacement.")),
function(x,a) stats::dhyper(x,a$success,a$failure,a$draws), function(x,a) stats::phyper(x,a$success,a$failure,a$draws), function(pr,a) stats::qhyper(pr,a$success,a$failure,a$draws), function(n,a) stats::rhyper(n,a$success,a$failure,a$draws),
function(a)c(max(0,a$draws-a$failure),min(a$draws,a$success)), function(a){N<-a$success+a$failure;m<-a$draws*a$success/N;v<-if(N>1)a$draws*(a$success/N)*(a$failure/N)*(N-a$draws)/(N-1) else 0;c(mean=m,variance=v,sd=sqrt(v))},
"The hypergeometric distribution is a classical finite-population model closely connected with combinatorial probability.",
"Used when sampling without replacement from a finite population containing a known number of successes and failures; it also underlies Fisher's exact test conditioning arguments.",
"Unlike the binomial distribution, draws are dependent because sampled items are not replaced.", NULL
)
reg$discreteuniform <- S(
"discreteuniform", "Discrete uniform", "discrete", "Discrete - bounded",
c("discrete uniform", "dunifint"), list(P("min", "Minimum integer", 1, -20, 20, 1, "Smallest possible integer."), P("max", "Maximum integer", 6, -20, 50, 1, "Largest possible integer.")),
function(x,a) ifelse(x==floor(x)&x>=a$min&x<=a$max,1/(a$max-a$min+1),0),
function(x,a) pmax(0,pmin(1,(floor(x)-a$min+1)/(a$max-a$min+1))),
function(pr,a) {
k <- a$max - a$min + 1
ifelse(pr <= 0, a$min, ifelse(pr >= 1, a$max,
pmax(a$min, pmin(a$max, a$min + ceiling(pr * k) - 1))))
},
function(n,a) sample(seq.int(a$min,a$max),n,replace=TRUE), function(a)c(a$min,a$max),
function(a){m<-(a$min+a$max)/2;k<-a$max-a$min+1;v<-(k^2-1)/12;c(mean=m,variance=v,sd=sqrt(v))},
"The discrete-uniform model formalizes equally likely values from a finite set, such as a fair die.",
"Useful for randomization examples, equally likely integer outcomes, and introductory probability teaching.",
"Every integer from min through max has exactly the same probability.", NULL
)
reg$betabinom <- S(
"betabinom", "Beta-binomial", "discrete", "Discrete - overdispersed binomial",
c("beta-binomial", "beta binomial"), list(P("n", "Number of trials", 10, 1, 200, 1, "Number of Bernoulli trials."), P("shape1", "Alpha", 2, .05, 20, .05, "First beta shape parameter."), P("shape2", "Beta", 2, .05, 20, .05, "Second beta shape parameter.")),
function(x,a){ok<-x==floor(x)&x>=0&x<=a$n;out<-rep(0,length(x));out[ok]<-exp(lchoose(a$n,x[ok])+lbeta(x[ok]+a$shape1,a$n-x[ok]+a$shape2)-lbeta(a$shape1,a$shape2));out},
function(x,a){vapply(x,function(xx){k<-0:min(a$n,floor(xx));if(xx<0)0 else if(xx>=a$n)1 else sum(exp(lchoose(a$n,k)+lbeta(k+a$shape1,a$n-k+a$shape2)-lbeta(a$shape1,a$shape2)))},numeric(1))},
function(pr,a){vapply(pr,function(pp){
if (pp <= 0) return(0)
if (pp >= 1) return(a$n)
k<-0:a$n
cdf<-cumsum(exp(lchoose(a$n,k)+lbeta(k+a$shape1,a$n-k+a$shape2)-lbeta(a$shape1,a$shape2)))
hit <- which(cdf >= pp)[1L]
if (is.na(hit)) a$n else k[hit]
},numeric(1))},
function(nn,a){pp<-stats::rbeta(nn,a$shape1,a$shape2);stats::rbinom(nn,a$n,pp)}, function(a)c(0,a$n),
function(a){p<-a$shape1/(a$shape1+a$shape2);rho<-1/(a$shape1+a$shape2+1);m<-a$n*p;v<-a$n*p*(1-p)*(1+(a$n-1)*rho);c(mean=m,variance=v,sd=sqrt(v))},
"The beta-binomial combines a binomial count with beta-distributed heterogeneity in the event probability.",
"Useful for overdispersed binomial counts, clustered binary outcomes, and teaching how heterogeneity in p inflates variance beyond the simple binomial model.",
"The induced intraclass-correlation/overdispersion parameter is rho = 1/(alpha+beta+1).", NULL
)
reg$zip <- S(
"zip", "Zero-inflated Poisson", "discrete", "Discrete - zero-inflated counts",
c("zero-inflated poisson", "zero inflated poisson"), list(P("lambda", "Poisson mean", 3, .01, 30, .05, "Mean of the Poisson count component."), P("zi", "Structural-zero probability", .2, 0, .95, .01, "Probability of an extra structural zero.")),
function(x,a){base<-(1-a$zi)*stats::dpois(x,a$lambda);base[x==0]<-a$zi+(1-a$zi)*exp(-a$lambda);base},
function(x,a) ifelse(x<0,0,a$zi+(1-a$zi)*stats::ppois(x,a$lambda)),
function(pr,a){m0<-a$zi+(1-a$zi)*exp(-a$lambda);vapply(pr,function(pp){if(pp<=m0)0 else stats::qpois((pp-a$zi)/(1-a$zi),a$lambda)},numeric(1))},
function(n,a){z<-stats::rbinom(n,1,a$zi);ifelse(z==1,0,stats::rpois(n,a$lambda))}, function(a)c(0,Inf),
function(a){m<-(1-a$zi)*a$lambda;v<-(1-a$zi)*a$lambda*(1+a$zi*a$lambda);c(mean=m,variance=v,sd=sqrt(v))},
"Zero-inflated count models were developed to represent data generated by a mixture of a structural-zero state and an ordinary count process.",
"Useful when there is a defensible two-process mechanism producing more zeros than a Poisson model can explain, for example structural non-risk plus a count process among those at risk.",
"A large number of zeros alone does not prove zero inflation; negative-binomial or hurdle models may be more appropriate depending on the data-generating process.", NULL
)
reg$zipf <- S(
"zipf", "Zipf (finite)", "discrete", "Discrete - heavy tailed",
c("zipfian"), list(P("exponent", "Exponent (s)", 1.5, .1, 5, .05, "Controls the power-law decay."), P("max", "Maximum rank", 50, 2, 500, 1, "Finite upper rank used by this teaching implementation.")),
function(x,a){k<-seq_len(a$max);den<-sum(k^(-a$exponent));ifelse(x==floor(x)&x>=1&x<=a$max,x^(-a$exponent)/den,0)},
function(x,a){k<-seq_len(a$max);cdf<-cumsum(k^(-a$exponent))/sum(k^(-a$exponent));vapply(x,function(xx)if(xx<1)0 else if(xx>=a$max)1 else cdf[floor(xx)],numeric(1))},
function(pr,a){k<-seq_len(a$max);cdf<-cumsum(k^(-a$exponent))/sum(k^(-a$exponent));vapply(pr,function(pp)k[which(cdf>=pp)[1L]],numeric(1))},
function(n,a){k<-seq_len(a$max);sample(k,n,replace=TRUE,prob=k^(-a$exponent))}, function(a)c(1,a$max),
function(a){k<-seq_len(a$max);w<-k^(-a$exponent);w<-w/sum(w);m<-sum(k*w);v<-sum((k-m)^2*w);c(mean=m,variance=v,sd=sqrt(v))},
"George Kingsley Zipf popularized an inverse-rank frequency law in linguistics and related empirical systems in the 20th century.",
"Used for rank-frequency and heavy-tailed discrete phenomena. This R4VN implementation is deliberately finite so probabilities, moments, and simulation are well defined for teaching.",
"The finite maximum is part of the model here; changing it changes the normalizing constant and probabilities.", NULL
)
reg$logseries <- S(
"logseries", "Logarithmic series", "discrete", "Discrete - positive counts",
c("log-series", "logarithmic"), list(P("p", "Parameter p", .6, .001, .99, .01, "Shape parameter between 0 and 1.")),
function(x,a) ifelse(x==floor(x)&x>=1,-a$p^x/(x*log(1-a$p)),0),
function(x,a){vapply(x,function(xx){if(xx<1)return(0);k<-seq_len(max(1,floor(xx)));min(1,sum(-a$p^k/(k*log(1-a$p))))},numeric(1))},
function(pr,a){vapply(pr,function(pp){
if (pp <= 0) return(1)
if (pp >= 1) return(Inf)
s<-0;k<-0L
while(s<pp&&k<100000L){k<-k+1L;s<-s-a$p^k/(k*log(1-a$p))}
as.numeric(k)
},numeric(1))},
function(n,a){u<-stats::runif(n);vapply(u,function(pp){s<-0;k<-0L;while(s<pp&&k<100000L){k<-k+1L;s<-s-a$p^k/(k*log(1-a$p))};k},integer(1))},
function(a)c(1,Inf), function(a){m<- -a$p/((1-a$p)*log(1-a$p));v<- -a$p*(a$p+log(1-a$p))/((1-a$p)^2*log(1-a$p)^2);c(mean=m,variance=v,sd=sqrt(v))},
"The logarithmic-series distribution was developed in early-20th-century probability and became well known through R. A. Fisher's work on species abundance.",
"Used for positive, strongly right-skewed counts and historically for species-frequency models.",
"Support starts at 1, unlike Poisson and negative-binomial count distributions which include zero.", NULL
)
reg
}
.r4vn_dist_find <- function(distribution) {
if (missing(distribution) || is.null(distribution) || !length(distribution)) stop("Supply a distribution name, for example `\"normal\"`, `\"binomial\"`, or `\"poisson\"`.", call. = FALSE)
z <- tolower(trimws(as.character(distribution)[1L]))
reg <- .r4vn_dist_registry()
hit <- names(reg)[vapply(reg, function(s) z %in% tolower(s$aliases), logical(1))]
if (!length(hit)) {
labels <- vapply(reg, `[[`, character(1), "id")
stop("Unknown distribution `", distribution, "`. Available names include: ", paste(labels, collapse = ", "), ".", call. = FALSE)
}
reg[[hit[1L]]]
}
.r4vn_dist_validate <- function(spec, a) {
finite1 <- function(nm, positive = FALSE, nonnegative = FALSE, integer = FALSE) {
z <- a[[nm]]
if (!is.numeric(z) || length(z) != 1L || !is.finite(z)) stop("Parameter `", nm, "` must be one finite number.", call. = FALSE)
if (positive && z <= 0) stop("Parameter `", nm, "` must be > 0.", call. = FALSE)
if (nonnegative && z < 0) stop("Parameter `", nm, "` must be >= 0.", call. = FALSE)
if (integer && abs(z-round(z)) > sqrt(.Machine$double.eps)) stop("Parameter `", nm, "` must be an integer.", call. = FALSE)
}
id <- spec$id
positive <- intersect(names(a), c("sd","sdlog","rate","shape","shape1","shape2","scale","df","df1","df2","size","lambda","exponent"))
for (nm in positive) finite1(nm, positive = TRUE)
if ("ncp" %in% names(a)) { finite1("ncp"); if (id %in% c("chisq","f") && a$ncp < 0) stop("`ncp` must be >= 0 for this distribution.", call. = FALSE) }
if ("p" %in% names(a)) { finite1("p"); if (a$p < 0 || a$p > 1 || (id %in% c("geometric","nbinom","logseries") && a$p <= 0) || (id=="logseries" && a$p>=1)) stop("Parameter `p` is outside its valid range.", call. = FALSE) }
if ("zi" %in% names(a)) { finite1("zi"); if (a$zi < 0 || a$zi >= 1) stop("`zi` must satisfy 0 <= zi < 1.", call. = FALSE) }
if (id %in% c("binomial","betabinom")) { finite1("n", nonnegative = TRUE, integer = TRUE); a$n <- as.integer(round(a$n)) }
if (id == "hypergeom") {
for(nm in c("success","failure","draws")) finite1(nm, nonnegative = TRUE, integer = TRUE)
if (a$draws > a$success + a$failure) stop("`draws` cannot exceed the population size.", call. = FALSE)
}
if (id %in% c("uniform","discreteuniform")) {
finite1("min"); finite1("max"); if (a$max <= a$min) stop("`max` must be greater than `min`.", call. = FALSE)
if (id=="discreteuniform" && (a$min!=round(a$min)||a$max!=round(a$max))) stop("Discrete-uniform bounds must be integers.", call. = FALSE)
}
if (id == "triangular") {
finite1("min"); finite1("mode"); finite1("max")
if (!(a$min < a$mode && a$mode < a$max))
stop("Require min < mode < max for the triangular distribution.", call. = FALSE)
}
if (id == "zipf") { finite1("max", positive=TRUE, integer=TRUE); a$max<-as.integer(round(a$max)) }
a
}
.r4vn_dist_args <- function(spec, dots) {
defaults <- setNames(lapply(spec$params, `[[`, "default"), vapply(spec$params, `[[`, character(1), "name"))
if (is.null(names(dots))) names(dots) <- rep("", length(dots))
unnamed <- which(!nzchar(names(dots)))
if (length(unnamed)) stop("Distribution parameters in `...` must be named.", call. = FALSE)
# Friendly aliases are deliberately distribution-specific so that a typo or
# an irrelevant parameter is never silently ignored.
allowed <- names(defaults)
move_alias <- function(old, new) {
if (old %in% names(dots)) {
if (new %in% names(dots))
stop("Supply only one of `", old, "` and `", new, "`.", call. = FALSE)
dots[[new]] <<- dots[[old]]
dots[[old]] <<- NULL
}
}
if ("p" %in% allowed) {
if ("prob" %in% names(dots)) move_alias("prob", "p")
if ("probability" %in% names(dots)) move_alias("probability", "p")
}
if (spec$id %in% c("binomial", "betabinom") && "trials" %in% names(dots))
move_alias("trials", "n")
if (spec$id == "binomial" && "size" %in% names(dots))
move_alias("size", "n")
if (spec$id == "gamma" && "scale" %in% names(dots)) {
if ("rate" %in% names(dots))
stop("For gamma, supply either `rate` or `scale`, not both.", call. = FALSE)
z <- dots$scale
if (!is.numeric(z) || length(z) != 1L || !is.finite(z) || z <= 0)
stop("Gamma `scale` must be one finite number > 0.", call. = FALSE)
dots$rate <- 1 / z
dots$scale <- NULL
}
if (spec$id == "nbinom" && "mu" %in% names(dots)) {
if ("p" %in% names(dots))
stop("For negative binomial, supply either `p` or `mu`, not both.", call. = FALSE)
mu <- dots$mu
size <- dots$size %||% defaults$size
if (!is.numeric(mu) || length(mu) != 1L || !is.finite(mu) || mu < 0)
stop("Negative-binomial `mu` must be one finite number >= 0.", call. = FALSE)
if (!is.numeric(size) || length(size) != 1L || !is.finite(size) || size <= 0)
stop("Negative-binomial `size` must be one finite number > 0.", call. = FALSE)
dots$p <- size / (size + mu)
dots$mu <- NULL
}
if (spec$id == "betabinom" && any(c("p", "rho") %in% names(dots))) {
if (!all(c("p", "rho") %in% names(dots)))
stop("For beta-binomial mean/correlation parameterization, supply both `p` and `rho`.", call. = FALSE)
if (any(c("shape1", "shape2") %in% names(dots)))
stop("For beta-binomial, use either `shape1`/`shape2` or `p`/`rho`, not both.", call. = FALSE)
pp <- dots$p; rho <- dots$rho
if (!is.numeric(pp) || length(pp) != 1L || !is.finite(pp) || pp <= 0 || pp >= 1 ||
!is.numeric(rho) || length(rho) != 1L || !is.finite(rho) || rho <= 0 || rho >= 1)
stop("For beta-binomial `p` and `rho`, require 0 < p < 1 and 0 < rho < 1.", call. = FALSE)
total <- 1 / rho - 1
dots$shape1 <- pp * total
dots$shape2 <- (1 - pp) * total
dots$p <- NULL
dots$rho <- NULL
}
extras <- setdiff(names(dots), allowed)
if (length(extras))
stop("Unused parameter(s) for ", spec$label, ": ", paste(extras, collapse = ", "), ".", call. = FALSE)
for (nm in intersect(names(dots), allowed)) defaults[[nm]] <- dots[[nm]]
.r4vn_dist_validate(spec, defaults)
}
.r4vn_dist_param_text <- function(spec, a, digits = 4L) {
paste(vapply(names(a), function(nm) paste0(nm, " = ", format(signif(a[[nm]], digits), trim=TRUE, scientific=FALSE)), character(1)), collapse = ", ")
}
.r4vn_dist_support_text <- function(spec, a) {
s <- spec$support(a)
lo <- if (is.infinite(s[1])) "-Inf" else format(s[1], trim=TRUE)
hi <- if (is.infinite(s[2])) "Inf" else format(s[2], trim=TRUE)
paste0(if(spec$type=="discrete") "integers in " else "", "[",lo,", ",hi,"]")
}
.r4vn_dist_probability_table <- function(spec, a, x, digits = 6L) {
if (is.null(x) || !length(x)) return(NULL)
x <- as.numeric(x)
rows <- list()
for (xx in x) {
if (spec$type == "discrete") {
eq <- spec$d(xx,a)
lt <- spec$p(ceiling(xx)-1,a)
le <- spec$p(floor(xx),a)
gt <- 1-le
ge <- 1-lt
} else {
eq <- 0
lt <- le <- spec$p(xx,a)
gt <- ge <- 1-le
}
rows[[length(rows)+1L]] <- data.frame(
x = xx,
Quantity = if(spec$type=="discrete") c("P(X = x)","P(X < x)","P(X <= x)","P(X > x)","P(X >= x)") else c("f(x) density","P(X < x)","P(X <= x)","P(X > x)","P(X >= x)"),
Value = c(if(spec$type=="discrete")eq else spec$d(xx,a),lt,le,gt,ge),
stringsAsFactors=FALSE
)
}
out <- do.call(rbind,rows)
# Keep full machine precision in the returned object. `digits` is a display
# control only; rounding here would make programmatic results differ from
# stats::dbinom()/pbinom() and related base-R reference functions.
out$Value <- pmax(0, pmin(ifelse(grepl("density", out$Quantity), Inf, 1), out$Value))
out
}
.r4vn_dist_range_table <- function(spec,a,lower,upper,digits=6L) {
if (is.null(lower) || is.null(upper)) return(NULL)
if (length(lower)!=1L||length(upper)!=1L||!is.finite(lower)||!is.finite(upper)||lower>upper) stop("Supply finite `lower <= upper`.",call.=FALSE)
if(spec$type=="discrete") {
inside <- spec$p(floor(upper),a)-spec$p(ceiling(lower)-1,a)
} else {
inside <- spec$p(upper,a)-spec$p(lower,a)
}
inside <- max(0,min(1,inside))
data.frame(Interval=c("P(lower <= X <= upper)","P(X outside interval)"),
lower=lower,upper=upper,Value=c(inside,1-inside),stringsAsFactors=FALSE)
}
.r4vn_dist_quantile_table <- function(spec,a,probs,digits=6L) {
if (is.null(probs) || !length(probs)) return(NULL)
probs <- as.numeric(probs)
if(any(!is.finite(probs)|probs<0|probs>1))stop("`probs` must contain probabilities between 0 and 1.",call.=FALSE)
data.frame(Probability=probs,Quantile=spec$q(probs,a),stringsAsFactors=FALSE)
}
.r4vn_dist_display_table <- function(x, digits = 6L, columns = NULL) {
if (is.null(x)) return(NULL)
out <- x
if (is.null(columns)) columns <- names(out)[vapply(out, is.numeric, logical(1))]
columns <- intersect(columns, names(out))
for (nm in columns) {
z <- out[[nm]]
ok <- is.finite(z)
z[ok] <- signif(z[ok], digits)
out[[nm]] <- z
}
out
}
.r4vn_dist_grid <- function(spec,a,x=NULL,lower=NULL,upper=NULL,points=501L) {
sup <- spec$support(a)
probs <- if(spec$type=="discrete") c(.001,.999) else c(.005,.995)
qq <- suppressWarnings(spec$q(probs,a))
lo <- if(is.finite(sup[1]))sup[1] else qq[1]
hi <- if(is.finite(sup[2]))sup[2] else qq[2]
extras <- c(x,lower,upper);extras<-extras[is.finite(extras)]
if(length(extras)){lo<-min(lo,extras);hi<-max(hi,extras)}
if(!is.finite(lo)||!is.finite(hi)||lo==hi){lo<- -5;hi<-5}
pad <- .05*(hi-lo); if(spec$type=="continuous"){lo<-if(is.finite(sup[1]))max(sup[1],lo-pad) else lo-pad;hi<-if(is.finite(sup[2]))min(sup[2],hi+pad) else hi+pad}
if(spec$type=="discrete") {
lo <- ceiling(lo); hi <- floor(hi)
if (is.finite(sup[1])) lo <- max(lo, ceiling(sup[1]))
if (is.finite(sup[2])) hi <- min(hi, floor(sup[2]))
if (!is.finite(lo)) lo <- 0
if (!is.finite(hi)) hi <- lo + 50
if (hi < lo) hi <- lo
if(hi-lo>250){center<-suppressWarnings(spec$q(.5,a));if(!is.finite(center))center<-(lo+hi)/2;lo<-max(lo,floor(center-125));hi<-min(hi,ceiling(center+125))}
seq.int(lo,hi)
} else seq(lo,hi,length.out=max(101L,as.integer(points)))
}
#' Explore a probability distribution and calculate probabilities
#'
#' `distdata()` is the command-line probability-distribution calculator for
#' R4VN. It combines distribution properties, point probabilities/densities,
#' tail probabilities, interval probabilities, quantiles, simulation, and a
#' Viewer-ready plot. [distlearn()] provides the interactive Shiny companion.
#'
#' @param distribution Distribution name. Common abbreviations are accepted,
#' for example `"normal"`, `"binomial"`, `"poisson"`, `"t"`, `"chisq"`,
#' `"weibull"`, `"betabinom"`, and `"zip"`.
#' @param ... Named parameters of the selected distribution. For example,
#' `n` and `p` for binomial, `mean` and `sd` for normal, or `lambda` for
#' Poisson. Defaults are teaching-friendly and are shown by `distlearn()`.
#' @param x Optional value(s). For discrete distributions R4VN reports
#' `P(X=x)`, `<`, `<=`, `>`, and `>=`. For continuous distributions it
#' reports density `f(x)` and the corresponding tail probabilities.
#' @param lower,upper Optional interval bounds for an interval probability.
#' @param probs Optional cumulative probabilities for which quantiles are
#' requested, for example `probs = c(.025, .5, .975)`.
#' @param nsim Optional number of random observations to simulate. `0` (the
#' default) does not simulate.
#' @param seed Optional user-supplied random seed for simulation. The default
#' `NULL` does not set a seed.
#' @param plot Logical; retain plot data and display the probability function
#' in the R4VN Viewer.
#' @param digits Number of significant digits in numerical probability output.
#' @param show Logical; open the R4VN Viewer result.
#' @param console Logical; also print the tabular result to the console.
#'
#' @details
#' A particularly useful teaching call is
#' `distdata("binomial", n = 10, x = 3, p = .2)`. With only these inputs R4VN
#' reports `P(X=3)`, `P(X<3)`, `P(X<=3)`, `P(X>3)`, and `P(X>=3)`, together
#' with the distribution's mean, variance, standard deviation, parameters,
#' support, and plot. Supplying `lower` and `upper` adds the probability inside
#' and outside an interval; supplying `probs` adds quantiles.
#'
#' The beta-binomial accepts either `shape1`/`shape2` or the more interpretable
#' pair `p`/`rho`. The negative binomial accepts `size` with either `p` or
#' `mu`. Gamma accepts `rate` or `scale`.
#'
#' @return Invisibly returns an object of classes `r4vn_distdata` and
#' `r4vn_stat`. Its `raw` component contains the distribution specification,
#' parameters, point/range probabilities, quantiles, simulation, and plot
#' grid.
#' @export
#'
#' @examples
#' distdata("binomial", n = 10, x = 3, p = .2, show = FALSE)
#' distdata("binomial", n = 20, p = .35, lower = 5, upper = 10, show = FALSE)
#' distdata("normal", mean = 100, sd = 15, x = 130, show = FALSE)
#' distdata("normal", mean = 100, sd = 15, probs = c(.025, .5, .975), show = FALSE)
#' distdata("poisson", lambda = 2.5, x = 0:4, show = FALSE)
#' distdata("nbinom", size = 2, mu = 5, x = 0:4, show = FALSE)
#' distdata("betabinom", n = 20, p = .3, rho = .1, x = 0:5, show = FALSE)
distdata <- function(distribution, ..., x = NULL, lower = NULL, upper = NULL,
probs = NULL, nsim = 0L, seed = NULL, plot = TRUE,
digits = 6L, show = TRUE, console = FALSE) {
call <- match.call()
if (!is.numeric(digits) || length(digits) != 1L || !is.finite(digits) || digits < 1 || digits > 15)
stop("`digits` must be one integer from 1 to 15.", call. = FALSE)
digits <- as.integer(round(digits))
spec <- .r4vn_dist_find(distribution)
dots <- list(...)
a <- .r4vn_dist_args(spec,dots)
moments <- spec$moments(a)
support <- .r4vn_dist_support_text(spec,a)
params_text <- .r4vn_dist_param_text(spec,a)
prop <- data.frame(
Property = c("Distribution","Type","Parameters","Support","Mean","Variance","Standard deviation"),
Value = c(spec$label,tools::toTitleCase(spec$type),params_text,support,
ifelse(is.na(moments["mean"]),"Not defined / not summarized",format(signif(moments["mean"],digits),trim=TRUE)),
ifelse(is.na(moments["variance"]),"Not defined / not summarized",format(signif(moments["variance"],digits),trim=TRUE)),
ifelse(is.na(moments["sd"]),"Not defined / not summarized",format(signif(moments["sd"],digits),trim=TRUE))),
stringsAsFactors=FALSE
)
point <- .r4vn_dist_probability_table(spec,a,x,digits)
rng <- .r4vn_dist_range_table(spec,a,lower,upper,digits)
quant <- .r4vn_dist_quantile_table(spec,a,probs,digits)
sections <- list("Distribution"=prop)
if(!is.null(point)) sections[["Probability at x"]] <- .r4vn_dist_display_table(point, digits, "Value")
if(!is.null(rng)) sections[["Interval probability"]] <- .r4vn_dist_display_table(rng, digits, "Value")
if(!is.null(quant)) sections[["Quantiles"]] <- .r4vn_dist_display_table(quant, digits, "Quantile")
sim <- NULL
if (!is.numeric(nsim) || length(nsim) != 1L || !is.finite(nsim) || nsim < 0 || abs(nsim - round(nsim)) > sqrt(.Machine$double.eps))
stop("`nsim` must be a non-negative integer.", call. = FALSE)
nsim <- as.integer(round(nsim))
if(nsim>0L){
if(!is.null(seed)){
seed <- .r4vn_seed_normalize(seed)
had<-exists(".Random.seed",envir=.GlobalEnv,inherits=FALSE);if(had)old<-get(".Random.seed",envir=.GlobalEnv,inherits=FALSE)
on.exit({if(had)assign(".Random.seed",old,envir=.GlobalEnv) else if(exists(".Random.seed",envir=.GlobalEnv,inherits=FALSE))rm(".Random.seed",envir=.GlobalEnv)},add=TRUE)
set.seed(seed)
}
sim<-spec$r(nsim,a)
sections[["Simulation"]]<-data.frame(n=nsim,Mean=signif(mean(sim),digits),SD=signif(stats::sd(sim),digits),Minimum=signif(min(sim),digits),Median=signif(stats::median(sim),digits),Maximum=signif(max(sim),digits),stringsAsFactors=FALSE)
}
grid <- if(isTRUE(plot)) .r4vn_dist_grid(spec,a,x,lower,upper) else NULL
raw <- list(spec=spec,parameters=a,moments=moments,point=point,range=rng,quantiles=quant,simulation=sim,grid=grid,x=x,lower=lower,upper=upper)
notes <- c(spec$use, spec$notes)
out <- .r4vn_result(paste0(spec$label," distribution"),sections,notes=notes,raw=raw,call=call)
class(out)<-c("r4vn_distdata",class(out))
.r4vn_show(out,show=show,console=console)
}
#' Plot an R4VN probability-distribution result
#'
#' @param x An object returned by [distdata()].
#' @param type `"density"`/`"mass"` or `"cdf"`.
#' @param main Optional plot title.
#' @param ... Additional graphical parameters passed to base graphics where
#' applicable.
#' @return The input object, invisibly.
#' @export
plot.r4vn_distdata <- function(x, type = c("density","cdf"), main = NULL, ...) {
type <- match.arg(type)
spec<-x$raw$spec;a<-x$raw$parameters;xx<-x$raw$grid
if(is.null(xx))xx<-.r4vn_dist_grid(spec,a,x$raw$x,x$raw$lower,x$raw$upper)
ttl<-main %||% paste0(spec$label,if(type=="cdf")" - cumulative distribution" else if(spec$type=="discrete")" - probability mass" else " - density")
if(type=="cdf"){
yy<-spec$p(xx,a)
if(spec$type=="discrete")graphics::plot(xx,yy,type="s",xlab="x",ylab="F(x)",main=ttl,...) else graphics::plot(xx,yy,type="l",xlab="x",ylab="F(x)",main=ttl,...)
graphics::abline(h=c(0,1),lty=3)
} else if(spec$type=="discrete"){
yy<-spec$d(xx,a);graphics::plot(xx,yy,type="h",lwd=3,xlab="x",ylab="P(X = x)",main=ttl,...);graphics::points(xx,yy,pch=16)
} else {
yy<-spec$d(xx,a);graphics::plot(xx,yy,type="l",lwd=2,xlab="x",ylab="Density",main=ttl,...)
}
if(!is.null(x$raw$x)){
xv<-x$raw$x[is.finite(x$raw$x)];if(length(xv))graphics::abline(v=xv,lty=2)
}
if(!is.null(x$raw$lower)&&!is.null(x$raw$upper))graphics::abline(v=c(x$raw$lower,x$raw$upper),lty=2)
invisible(x)
}
.r4vn_dist_svg <- function(x,type="density") {
f<-tempfile(fileext=".svg")
grDevices::svg(f,width=7.5,height=4.3,pointsize=11)
on.exit({try(grDevices::dev.off(),silent=TRUE);unlink(f)},add=TRUE)
plot(x,type=type)
grDevices::dev.off()
txt<-readLines(f,warn=FALSE,encoding="UTF-8")
start<-grep("<svg",txt,fixed=TRUE)[1L]
if(is.na(start))return("")
paste(txt[start:length(txt)],collapse="\n")
}
.r4vn_viewer_distdata <- function(x) {
spec<-x$raw$spec
hero<-paste0('<section class="r4vn-section"><div style="display:grid;grid-template-columns:1fr 1fr;gap:12px">',
'<div><h2>Why this distribution?</h2><p>',.r4vn_view_escape(spec$use),'</p></div>',
'<div><h2>Historical note</h2><p>',.r4vn_view_escape(spec$history),'</p></div></div></section>')
graph<-if(is.null(x$raw$grid))"" else paste0('<section class="r4vn-section"><h2>Distribution plot</h2>',.r4vn_dist_svg(x,"density"),'</section>')
blocks<-paste0(vapply(names(x$sections),function(nm).r4vn_view_section(nm,x$sections[[nm]]),character(1)),collapse="")
formula<-if(is.null(spec$formula)||!nzchar(spec$formula))"" else paste0('<section class="r4vn-section"><h2>Formula / definition</h2><pre>',.r4vn_view_escape(spec$formula),'</pre></section>')
.r4vn_view_document(x$title,paste0(hero,graph,blocks,formula),notes=x$notes,subtitle="R4VN distribution calculator",prefix="r4vn-distdata-")
}
Try the R4VN 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.