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
______________________________________________
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.
