[R] The Percentile of a User-Defined pdf

[Nissim Kaufmann]

I would like to give a probability distribution function of a function of
(x,y) on the half-plane y>0, and a constant 0<c<1 and I would like to know
the c percentile of the marginal distribution of x.  I have tried along
the lines of the following but I keep getting errors:

# SIMPLIFIED PROBLEM
# The plan is to solve for the .975 percentile "xc" of the marginal x
distribution of the pdf (say it is proportional to 1/(1+x^2+y^2) for
simplicity) which has support on the real plane.
# The function 1/(1+x^2+y^2) has value (normalization constant)
approximately equal to "I" that I was able to
# program with no problem, as shown below.

# Approximate I, the normalization constant.
# This works fine.
c=1e+3 #the bound of integration in the three directions; if I use Inf I
get error.

I=integrate(function(y) {
   sapply(y, function(y) {
     integrate(function(x) {
       sapply(x, function(x) 1/(1+x^2+y^2))
     }, -c, c)$value
   })
 }, -c, c)

# Preliminary step -- define function J as an integral up to variable xc.
# I am still stuck on this step -- R says
# "Error in is.vector(X): object 'y' not found."

J=sapply(xc, function(xc) {integrate(function(x) {
   sapply(y, function(x) {
     integrate(function(y) {
       sapply(x, function(y) 1/(1+x^2+y^2))
     }, -c, c)$value
   })
  }, -c, xc)$value
 })

# Final step -- solve for .975 percentile of the above function J
uniroot(
sapply(xc, function(xc) {integrate(function(x) {
   sapply(y, function(x) {
     integrate(function(y) {
       sapply(x, function(y) 1/(1+x^2+y^2))
     }, -c, c)$value
   })
  }, -c, xc)$value
 })
 )/I-.975,
lower=-c, upper=c, tol=1e-10)$root

I don't have much programming experience.  Thank you for your help.

Nissim Kaufmann
University at Albany