Tools for working with the Voigt distribution

Description

Random generation, density function and parameter estimation for the Voigt distribution. The main objective of this package is to provide R users with efficient estimation of Voigt parameters using classic iid data in a Bayesian framework. The estimating function allows flexible prior specification, specification of fixed parameters and several options for MCMC posterior simulation. A basic version of the algorithm is described in: Cannas M. and Piras, N. (2025) <doi:10.1007/978-3-031-96303-2_53>.

Details

The package contains functions for working with the Voigt distribution in the R environment:

rvoigt generates random variates

dvoigt calculates voigt density

evoigt estimates parameters following a Bayesian approach

otherwise unavailable within R (R < 4.4.1). The main function evoigt provides both point and interval estimates for the Voigt parameters using a Bayesian approach.

Author(s)

Massimo Cannas [aut, cre], Nicola Piras [aut]

Maintainer: Massimo Cannas <massimo.cannas@unica.it>

References

Kendall, D. G. (1938). The effect of radiation damping and Doppler broadening on the atomic absorption coefficient. Zeitschrift fur Astrophysik, Vol. 16, p.308 https://adsabs.harvard.edu/full/1938ZA.....16..308K

Cannas, M. and Piras, N. Mixture representation and parameter estimation for the Voigt profile (submitted)

Examples

## Not run: 
##  See examples in the help file of each function.
    
## End(Not run)

The Voigt Distribution

Description

Density function and random generation for the Voigt distribution with median equal to mu and scale parameters equal to sigma and gamma.

Usage

dvoigt(y, mu = 0, sigma = 1, gamma = 1, log = FALSE, complex = FALSE)
rvoigt(n = 1, mu = 0, sigma = 1, gamma = 1)

Arguments

Details

The Voigt distribution is the convolution of a Normal and a Cauchy. The density function is

f(y)=\frac{\Re(w(z))}{\sigma \sqrt{2\pi}}, \quad \sigma,\gamma>0

where w(z)=e^{-z^2}erfc(-iz) is the Faddeeva function, z=\frac{y-\mu +i\gamma}{\sigma\sqrt{2}}, mu is the median, sigma is the standard deviation of the normal component and gamma is the scale parameter of the Cauchy component. If mu, sigma and gamma are not specified they assume the default values of 0, 1 and 1, respectively.

Value

dvoigt gives the density and rvoigt generates random deviates.

Note

rvoigt is a wrapper of both rnorm and rcauchy functions. dvoigt uses erfz, the complex error function, see erfz.

Author(s)

Massimo Cannas [aut, cre], Nicola Piras [aut]

References

Kendall, D. G. (1938). The effect of radiation damping and Doppler broadening on the atomic absorption coefficient. Zeitschrift fur Astrophysik, Vol. 16, p.308 https://adsabs.harvard.edu/full/1938ZA.....16..308K

Abramowitz, M., and I. A. Stegun (1972). Handbook of Mathematical Functions (with Formulas, Graphs, and Mathematical Tables). Dover, New York. https://personal.math.ubc.ca/~cbm/aands/notes.htm

See Also

evoigt

Examples

dvoigt(0)
# 0.2087093
# plot voigt and gaussian density
x = seq(from=-4,to=4, by=0.01)
plot(x, dvoigt(x), type="l",ylab="",ylim=c(0,0.8))
lines(x,dvoigt(x,0,1/3,1/3),col="blue")
lines(x,dnorm(x),lty=2, col="red")
mtext("dvoigt(x,0,1,1)", adj = 0,col="black")
mtext("dnorm(x)", adj = 1/2,col="red")
mtext("dvoigt(x,0,1/3,1/3)", adj = 1,col="blue")

# compare true and empirical density
rvoigt(1)
x = rvoigt(n=500)
q=quantile(x,0.99)
hist(x[x>-q & x<q], prob=TRUE, breaks=30,ylim=c(0,1.1*dvoigt(0)),main="", xlab="x")
x.grid = seq(-q,q,by=.1)
lines(x.grid, dvoigt(x.grid), type="l", ylab="", xlab="",col="red")
# compare with cauchy and normal density
par(mfrow=c(1,2))
x = rvoigt(n=500, mu = 0,sigma = 1,gamma = 0.1)
q=quantile(x,0.99)
hist(x[x>-q & x<q], prob=TRUE,breaks=20,col="lightgreen",main = "dvoigt(0, 1, 0.1)",xlab="x")
curve(dnorm(x),lty=1, add=TRUE)

y = rvoigt(n=500, mu = 0,sigma = 0.1,gamma = 1)
q=quantile(y,0.99)
hist(y[y>-q & y<q], prob=TRUE,breaks=25,col="lightblue",main = "dvoigt(0, 0.1, 1)",xlab="x")
curve(dcauchy(x),lty=1, add=TRUE)
dev.off()

Estimation of Voigt distribution parameters

Description

The function provides both point and interval estimates for the Voigt distribution parameters in a Bayesian way (see Details).

Usage

evoigt(data, hyper = NULL, init = NULL, S = 10000, burn = S/2,
thin = 10, fix.arg = NULL, chain = FALSE, ...)

Arguments

Details

The function runs a Gibbs sampler for estimating parameters using the following prior distributions: \mu\sim N(\mu_0, \nu_0^2), \sigma\sim \Gamma(a, b) and \gamma\sim \Gamma(c, d).

Value

evoigt returns a list with the following components:

Author(s)

Massimo Cannas [aut, cre], Nicola Piras [aut]

References

Cannas, M. and Piras, N. (2025) Estimation of Voigt Distribution Parameters: A Bayesian Approach. (conference paper) https://link.springer.com/chapter/10.1007/978-3-031-96303-2_53

Cannas, M. and Piras, N. Mixture representation and parameter estimation for the Voigt profile (submitted)

See Also

See HPDinterval for details on how the highest posterior density interval is built.

Examples

x = rvoigt(500, mu=0, sigma=1,gamma=1)
# point estimates and (default) 95% credibility intervals
evoigt(data=x)
# point estimates and 90% credibility intervals
evoigt(data=x, prob = 0.9)
# chain values
res = evoigt(data=x, prob = 0.9, chain=TRUE)
mu.chain = res$mu.chain
plot(0: (length(mu.chain)-1), type = "l", mu.chain, xlab="",ylab=c(expression(mu)) )
points(0,mu.chain[1],pch=1,col="red")
# if a parameter is known its value can be fixed using "fix.arg"
# e.g. set sigma =1:
res = evoigt(data=x, fix.arg=c(NA, 1, NA))
res["posterior mean"]
res["posterior median"]
res["HPD interval"]