Skip to content

The Percentile of a User-Defined pdf

3 messages · Nissim Kaufmann, Dieter Menne, David Winsemius

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
Dieter Menne
#
Nissim Kaufmann wrote:
Once you are inside the first "{", R only knows about the x it received in
as a parameter (assuming there is no global y). I have not checked the
details of the code (it looks overly complicated), but 

   sapply(x, function(y) {

could work. In theory, keeping function(x) would also work, but I suggest
that you use another variable name in inner nests to avoid confusion (human,
not CPU).

My favorite step in developing nested xapply is to first remove all code and
only put a 

print(str(x))

inside the function, so I can clearly see what is passed in.

Dieter
David Winsemius
#
On Jan 1, 2011, at 10:42 PM, Nissim Kaufmann wrote:

            

      
      
      
    

        
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.

            

      
      
      
    

        
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.

            

      
      
      
    

        
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

            

      
      
      
    

        

      
      
    

        
David Winsemius, MD
West Hartford, CT