Reset function R
I am creating a set of functions that return probability density functions (pdfs) from PDFs of two independent random variables.
The most common example of this can be the sum of independent random variables X , Y , which is given by the convolution of their pdf files.
After this post, I defined the following function, which takes a couple of PDFs as an argument and returns them as a fold:
dSumXY <- function(dX, dY){
# Create convolution of distributions.
dReturn <- function(z){
integrate( function(x,z){
dX(x) * dY(z - x) },
-Inf, Inf, z)$value
}
# Vectorize convolution.
dReturn <- Vectorize(dReturn)
return(dReturn)
}
This works as expected in the following example:
# Define pdfs of two (identical) uniform [-1,1] distributions
unifX <- function(x) dunif(x, min = -1, max = 1)
unifY <- function(x) dunif(x, min = -1, max = 1)
# Find convolution of their pdfs.
convXY <- dSumXY(unifX, unifY)
# Plot the convolved pdf.
plot(seq(-3,3,by = 0.1), convXY(seq(-3,3,by = 0.1) ), type = 'l')
# Sample from the distribution
convSample <- runif(10000, min = -1, max = 1) + runif(10000, min = -1, max = 1)
# Plot density of sample.
lines( density(convSample) , col = "red" )
More generally, this works for a number of combinations of pairs of uniform distributions, however, when I try to overthrow a pair of uniform distributions [1,2], I don't get the true result:
# Define pdfs of two (identical) uniform [1,2] distributions
unifX2 <- function(x) dunif(x, min = 1, max = 2)
unifY2 <- function(x) dunif(x, min = 1, max = 2)
# Find convolution of their pdfs.
convXY2 <- dSumXY(unifX2, unifY2)
# Plot the convolved pdf.
plot(seq(1,5,by = 0.1), convXY2(seq(1,5,by = 0.1) ), type = 'l')
# Sample from the distribution
convSample2 <- runif(10000, min = 1, max = 2) + runif(10000, min = 1, max = 2)
# Plot density of sample.
lines( density(convSample2) , col = "red" )
In particular (with a bit of probabilistic knowledge!) It is clear that the pdf of the sum of two uniform [1,2] variables must be nonzero by 3.75; in particular, it must be 1/4. However
# This should be equal to 1/4, but returns 0.
convXY2(3.75)
I've tried this on two separate machines and copied the same problem, so I'm curious to know where the problem came from.
Thanks in advance.
source to share
dunif deals with integers ...
dunif(.9, 1, 2)
[1] 0
dunif(1, 1, 2)
[1] 1
dunif(1.5, 1, 2)
[1] 1
dunif(2, 1, 2)
[1] 1
dunif(2.1, 1, 2)
[1] 0
If you only use integers for your function and then interpolate, you will get the expected response.
convXY2(2)
[1] 0
convXY2(3)
[1] 1
convXY2(4)
[1] 0
z <- 3.75
remainder <- z %% 1
convXY2(floor(z))*(1-remainder) + convXY2(ceiling(z))*(remainder)
[1] 0.25
Your function works correctly for other contiguous distributions, so maybe you just need to account for that particular situation in the wrapper function. This is or use something other than dunif that effectively divides your probability based on the precision you give it.
dunif2 <- function(x, min, max) { if (x >= min & x <= max) { return(10^-nchar(strsplit(as.character(x), "\\.")[[1]][[2]])) } else { return(0) } }
dunif2(1.45, 1, 2)
[1] 0.01
dunif2(1.4, 1, 2)
[1] 0.1
source to share