[R] The Percentile of a User-Defined pdf

[David Winsemius]

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