The Percentile of a User-Defined pdf

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) {
It's a bad idea to define a new function "I". There already is one in  
R. It might or might not cause problems in this case, depending on  
whether any of the internal functions might be calling it.
  sapply(y, function(y) {
    integrate(function(x) {
      sapply(x, function(x) 1/(1+x^2+y^2))
    }, -c, c)$value
It's also a bad idea to use "-c" and "c" as your name for limits of  
integration (or for any other variable). "c" is a rather fundamental R  
function. It may not confuse the interpreter but it will at the vary  
least confuse the humans who attempt to understand it and will give  
error messages that are difficult to interpret.
  })
}, -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
})
In addition to the debugging strategy suggested by Menne, from the  
console you can issue a traceback() call immediately after a call and  
see the sequence of calls up to the error. You can also place a  
browser() call inside the sapply calls which will give you the  
capability of inspecting the local environment inside the call.

?browser
# 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

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
David Winsemius, MD
West Hartford, CT

The Percentile of a User-Defined pdf

Thread (3 messages)