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@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
______________________________________________
R-help@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.
