How to create a distribution function in R?

问题:

Given the following function:

f(x) = (1/2*pi)(1/(1+x^2/4))

How do I identify it's distribution and write this distribution function in R?

回答1:

So this is your function right now (hopefully you know how to write an R function; if not, check writing your own function):

f <- function (x) (pi / 2) * (1 / (1 + 0.25 * x ^ 2))

f is defined on (-Inf, Inf) so integration on this range gives an indefinite integral. Fortunately, it approaches to Inf at the speed of x ^ (-2), so the integral is well defined, and can be computed:

C <- integrate(f, -Inf, Inf)
# 9.869604 with absolute error < 1e-09

C <- C$value  ## extract integral value
# [1] 9.869604

Then you want to normalize f, as we know that a probability density should integrate to 1:

f <- function (x) (pi / 2) * (1 / (1 + 0.25 * x ^ 2)) / C

You can draw its density by:

curve(f, from = -10, to = 10)

回答2:

Now that I have the probably distribution function I was wondering how to create a random sample of say n = 1000 using this new distribution function?

An off-topic question, but OK to answer without your making a new thread. Useful as it turns out subtle.

Compare

set.seed(0); range(simf(1000, 1e-2))
#[1] -56.37246  63.21080
set.seed(0); range(simf(1000, 1e-3))
#[1] -275.3465  595.3771
set.seed(0); range(simf(1000, 1e-4))
#[1] -450.0979 3758.2528
set.seed(0); range(simf(1000, 1e-5))
#[1] -480.5991 8017.3802

So I think e = 1e-2 is reasonable. We could draw samples, make a (scaled) histogram and overlay density curve:

set.seed(0); x <- simf(1000)
hist(x, prob = TRUE, breaks = 50, ylim = c(0, 0.16))
curve(f, add = TRUE, col = 2, lwd = 2, n = 201)