How to find pdf of a joint distribution in R? $F(x,y) =\frac{1}{6}(x^2\, y+x\, y^2)\,,\quad 0\leq x\leq 2,\, 0\leq y\leq 1$
Above is the joint distribution given,


*

*how to find out cumulative distribution function of y?

*how to obtain joint probability density function of x and y?


I am a beginner in R, I know basic commands.
thanks for the help
 A: This is a double integral in R (It's not done symbolically as Mathematica would do it but rather numerically):
 llimy <- 0;  llimx=0
 ulimy <-  2    ; ulimx=1

 f <- function(x,y) 1/6*(x^2*y+x*y^2)


 integrate(function(y) { 
    sapply(y, function(y) {
      integrate(function(x) f(x,y), llimx, ulimx)$value
    })
  }, llimy, ulimy)
#  0.3333333 with absolute error < 3.7e-15

The joint pdf is just the function, f, divided by the value of the integral over the full range of the values.
A: If you're trying to do this symbolically, you may want to try Wolfram Alpha. If you don't understand how to do this symbolically by hand, neither Wolfram Alpha nor this post will help you. You'll need to consult your statistics textbook for that. What this post does answer is how to get R to numerically compute distributions.
The volume under the curve is 1/3, so we just multiply by 3 to get the probability distribution for x and y. Obviously R doesn't deal with symbolic algebra (without the Ryacas package), but it is fairly easy to make pdfs and cdfs of functions. There is probably a simpler or more computationally efficient way, but this solution is fast enough for what you may be trying to do.
First, we input the pdf of x and y.
pdfxy <- function(x, y) (x^2 * y + x * y^2)/2

We convert this to a pdf of just y by integrating over the possible x values. The sapply function makes it so this function can easily take vectors as the y argument.
pdfy <- function(y) {
  result <- sapply(y, function(b) integrate(function(a) pdfxy(a, b), 0, 2)$value)
  return(result)
}

We then make this into a cdf by integrating over y from 0 to the desired value. The sapply function in this function is only required if you want to be able to input vectors.
cdfy <- function(y) {
  result <- sapply(y, function(a) integrate(function(b) pdfy(b), 0, a)$value)
  return(result)
}

Now we can just type
> cdfy(c(0,0.5,1))
[1] 0.0000000 0.2083333 1.0000000
> cdfy(0.4)
[1] 0.128

You now have a function in R that calculates the cumulative distribution of y given any y.
