Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Estimating moments from a distribution function

I have a non-normal distribution function which I want to calculate it's moments (mean, variance, skewness, and kurtosis). The packages I know are e1071 and moments which calculate the moments for a vector of discrete values. Is there a package that estimates moments for a continuous distribution function? As an example assume I have this distribution:

tmp <- runif(100,-10,10)
PDFdata <- density(tmp , kernel = "gaussian")
PDF <- splinefun(PDFdata$x , PDFdata$y)

Now I want to calculate:

enter image description here

like image 384
Novic Avatar asked Mar 21 '26 16:03

Novic


2 Answers

You can use function integrate to compute the moments.
All you need is to define an integrand, which I call int in the code below.

First, I will make the results reproducible by setting the RNG seed.

set.seed(6126)    # Make the results reproducible

tmp <- runif(100,-10,10)
PDFdata <- density(tmp , kernel = "gaussian")
PDF <- splinefun(PDFdata$x , PDFdata$y)

int <- function(x, c = 0, g = PDF, n = 1) (x - c)^n * g(x)

integrate(int, lower = -15, upper = 15, n = 1)
#-0.3971095 with absolute error < 7.8e-06

integrate(int, lower = -15, upper = 15, n = 2)
#35.76295 with absolute error < 0.0012

plot(PDFdata)
curve(PDF, -15, 15, add = TRUE, col = "blue", lty = "dashed")

density plot

like image 102
Rui Barradas Avatar answered Mar 23 '26 05:03

Rui Barradas


If you are estimating your density from data, you're better off using empirical moments from the data to estimate the moments of the distribution. If you just used this as an example of a function, then you could use the integrate function from the stats package. For example,

set.seed(123)
tmp <- runif(100,-10,10)
PDFdata <- density(tmp , kernel = "gaussian")
PDF <- splinefun(PDFdata$x , PDFdata$y)

mean(tmp)
[1] -0.02882012
integrate(function(x) x*PDF(x), -Inf, Inf)
Error in integrate(function(x) x * PDF(x), -Inf, Inf) : 
  the integral is probably divergent

and indeed, for this dataset that PDF function is not a density:

plot(PDF, from = -100, to = 100)

enter image description here

So we force it to zero outside the range of the density estimate:

PDF2 <- function(x) ifelse(x < min(PDFdata$x) | x > max(PDFdata$x), 0, 
                           PDF(x))

integrate(function(x) x*PDF2(x), -Inf, Inf)
-0.02900585 with absolute error < 8.2e-05
like image 30
user2554330 Avatar answered Mar 23 '26 05:03

user2554330



Donate For Us

If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!