library(edfun) library(knitr) knitr::opts_chunk$set( # cache = TRUE, dpi = 60, comment = "#>", tidy = FALSE) # http://stackoverflow.com/questions/24091735/why-pandoc-does-not-retrieve-the-image-file # < ! -- rmarkdown v1 -->
Author: Tal Galili ( Tal.Galili@gmail.com )
As mentioned in CRAN Task View: Probability Distributions
Empirical distribution : Base R provides functions for univariate analysis: (1) the empirical density (see density()), (2) the empirical cumulative distribution function (see ecdf()), (3) the empirical quantile (see quantile()) and (4) random sampling (see sample()).
This package aims to easily wrap these into a single function edfun
(short for Empirical Distribution FUNctions). Also, since quantile is generally a slow function to perform, the default for creating a quantile function (inverse-CDF) is by approximating the function of predicting the data values (x) from their quantiles (CDF). This is done using the approxfun
function. It takes a bit longer to create qfun, but it is MUCH faster to run than quantile (and is thus much better for simulations). Special care is taken for dealing with the support of the distribution (if it is known upfront).
The added speed allows to use these functions to run simulation studies for unusual distributions.
To install the stable version on CRAN:
# install.packages('edfun') # not on CRAN yet
To install the GitHub version:
# You'll need devtools if (!require(devtools)) install.packages("devtools"); devtools::install_github('talgalili/edfun')
And then you may load the package using:
library("edfun")
First load the library:
library(edfun)
set.seed(123) x <- rnorm(1000) x_dist <- edfun(x) f <- x_dist$dfun curve(f, -2,2) f <- x_dist$pfun curve(f, -2,2) f <- x_dist$qfun curve(f, 0,1) new_x <- x_dist$rfun(1000) hist(new_x)
The same can be done if we had the density of the distribution. In this case, x
should be a sequence of values for which we wish to pre-compute the quantiles (inv-CDF).
x <- seq(-4,4, length.out = 1e3) x_dist <- edfun(x, dfun = dnorm) f <- x_dist$dfun curve(f, -2,2) f <- x_dist$pfun curve(f, -2,2) f <- x_dist$qfun curve(f, 0,1) new_x <- x_dist$rfun(1000) hist(new_x)
If dfun
is NULL
then it is just returned as NULL
. This is useful if using the functions ONLY for creating the CDF and inv-CDF and wanting to save computation time.
set.seed(123) x <- rnorm(1000) x_dist <- edfun(x, dfun = NULL) x_dist$dfun # but this still works just fine: f <- x_dist$pfun curve(f, -2,2) f <- x_dist$qfun curve(f, 0,1)
If it is known, we can specifically define the support:
x <- seq(-4,4, length.out = 1e3) x_dist_no_support <- edfun(x, dfun = dnorm) x_dist_known_support <- edfun(x, dfun = dnorm, support = c(-Inf, Inf)) x_dist_no_support$qfun(0) x_dist_known_support$qfun(0)
Timing is basically the same when using the functions relying on density (i.e.: dfun
) or on a sample only. But the accuracy of using the density is much better.
set.seed(2016-08-18) x_simu <- rnorm(1e5) x_funs_simu <- edfun(x_simu, support = c(-Inf, Inf)) x <- seq(-4,4, length.out = 1e5) x_funs <- edfun(x, dfun = dnorm, support = c(-Inf, Inf)) x_funs$qfun(0) # -Inf # precision - qfun q_to_test <- seq(0.001,.999, length.out = 100) edfun_out <- x_funs$qfun(q_to_test) # -4 edfun_simu_out <- x_funs_simu$qfun(q_to_test) # -4 real_out <- qnorm(q_to_test) mean(abs(edfun_out-real_out)) mean(abs(edfun_simu_out-real_out)) # quite terrible compared to when knowing dfun library(microbenchmark) microbenchmark(dfun_known = x_funs$qfun(q_to_test), dfun_NOT_known = x_funs_simu$qfun(q_to_test)) # same time for the q function! # precision - pfun q_to_test <- seq(-3,3, length.out = 100) edfun_out <- x_funs$pfun(q_to_test) # -4 edfun_simu_out <- x_funs_simu$pfun(q_to_test) # -4 real_out <- pnorm(q_to_test) mean(abs(edfun_out-real_out)) mean(abs(edfun_simu_out-real_out)) # quite terrible compared to when knowing dfun library(microbenchmark) microbenchmark(dfun_known = x_funs$pfun(q_to_test), dfun_NOT_known = x_funs_simu$pfun(q_to_test)) # same time for the p function! # timing for the rfun library(microbenchmark) microbenchmark(dfun_known = x_funs$rfun(100), dfun_NOT_known = x_funs_simu$rfun(100)) # this rfun is slower...
First load the library:
library(edfun)
# credit: library(smoothmest) ddoublex <- function (x, mu = 0, lambda = 1) { exp(-abs(x - mu)/lambda)/(2 * lambda) } curve(ddoublex, -4,4) # the peak in 0 should go to Inf, but it doesn't because of the limits of `curve` x <- seq(-4,4, length.out = 1e3) x_dist <- edfun(x, dfun = ddoublex) f <- x_dist$dfun curve(f, -4,4) f <- x_dist$pfun curve(f, -4,4) f <- x_dist$qfun curve(f, 0,1) new_x <- x_dist$rfun(1000) hist(new_x)
First load the library:
library(edfun)
# http://stackoverflow.com/questions/20452650/writing-a-bimodal-normal-distribution-function-in-r?rq=1 # random sample: # mixtools::rnormmix # random sample of a mixture of distributions # http://stats.stackexchange.com/questions/70855/generating-random-variables-from-a-mixture-of-normal-distributions # nor1mix::rnorMix # https://en.wikipedia.org/wiki/Mixture_distribution # density: dmixnorm <- function(x, p1 = 0.5, m1=0, m2=0, s1=1, s2=1) { p1 * dnorm(x, m1, s1) + (1-p1) * dnorm(x, m2, s2) } # a convex mixture of densities is a density: # https://en.wikipedia.org/wiki/Mixture_distribution#Properties # yay - it works. dmixnorm_1 <- function(x) dmixnorm(x, .75, 0,5, 1,1) curve(dmixnorm_1, -4,9) x <- seq(-4,9, length.out = 1e3) x_dist <- edfun(x, dfun = dmixnorm_1) f <- x_dist$dfun curve(f, -4,9) f <- x_dist$pfun curve(f, -4,9) f <- x_dist$qfun curve(f, 0,1) new_x <- x_dist$rfun(1000) hist(new_x)
You are welcome to:
You can see the most recent changes to the package in the NEWS.md file
sessionInfo()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.