R help - Jan 2011 - 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) {
   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

Nissim Kaufmann wrote:
> 
> 
> 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
>  })
> 
> 
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







            
-- 
View this message in context:
http://r.789695.n4.nabble.com/The-Percentile-of-a-User-Defined-pdf-tp3170672p3170786.html
Sent from the R help mailing list archive at Nabble.com.

On Jan 1, 2011, at 10:42 PM, Nissim Kaufmann wrote:
> 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

I got it to work:
# To get a percentile of a single-variable function:
# Step 1: Integrate over the domain to ge the normalization constant:
Z<-integrate(function(x) sqrt(1+x^-1), 1,2)$value
Z
# Step 2: Find the .975 percentile
x975<-uniroot(function(t) integrate(function(x) sqrt(1+x^-1), 1,
t)$value/Z-.975,
               lower=1, upper=2, tol=5e-4 )$root
x975

# To get a percentile of a marginal of a bivariate function of (x,y):
# Define and compute the marginal x distribution:
fx<-function(x) {sapply(x, function(x)
                                integrate( function(y)  sqrt(sin(x)+1/y),
0,10)$value
                                     )
                         }
# Then proceed as above in the singe-variable case.


Thank you Dieter and David.

Nissim Kaufmann
Dept. of Mathematics and Statistics
University at Albany

In reply to
http://www.mail-archive.com/r-help at r-project.org/msg121420.html