Path: news.mathworks.com!newsfeed-00.mathworks.com!newscon02.news.prodigy.net!prodigy.net!news.glorb.com!postnews.google.com!c4g2000hsg.googlegroups.com!not-for-mail
From: NZTideMan <mulgor@gmail.com>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Calculating cumulative probability
Date: Thu, 14 Feb 2008 14:50:14 -0800 (PST)
Organization: http://groups.google.com
Lines: 101
Message-ID: <53ab9842-75a9-4c46-82f2-f752beba4735@c4g2000hsg.googlegroups.com>
References: <fp1qjo$ka8$1@fred.mathworks.com> <fp24vh$49a$1@fred.mathworks.com> 
NNTP-Posting-Host: 202.78.152.105
Mime-Version: 1.0
Content-Type: text/plain; charset=ISO-8859-1
Content-Transfer-Encoding: quoted-printable
X-Trace: posting.google.com 1203029415 3543 127.0.0.1 (14 Feb 2008 22:50:15 GMT)
X-Complaints-To: groups-abuse@google.com
NNTP-Posting-Date: Thu, 14 Feb 2008 22:50:15 +0000 (UTC)
Complaints-To: groups-abuse@google.com
Injection-Info: c4g2000hsg.googlegroups.com; posting-host=202.78.152.105; 
User-Agent: G2/1.0
X-HTTP-UserAgent: Mozilla/4.0 (compatible; MSIE 7.0; Windows NT 5.1; Maxthon; 
Xref: news.mathworks.com comp.soft-sys.matlab:451509


On Feb 15, 10:53=A0am, "Roger Stafford"
<ellieandrogerxy...@mindspring.com.invalid> wrote:
> "Roger Stafford" <ellieandrogerxy...@mindspring.com.invalid> wrote in
> message <fp24vh$49...@fred.mathworks.com>...
>
>
>
> > =A0 Since you know X and Y to be independent random variables, your desi=
red
> > probability can be expressed as the iterated integral:
>
> > =A0P(k) =3D integral, 0 to inf, (integral, y to y+k, f(x)*g(y) dx) dy
>
> > where f(x) and g(y) are the respective probability density functions of =
X and
> Y. =A0
> > That is: the integral from 0 to infinity of the integral from y to y+k o=
f the
> > product f(x)*g(y), first with respect to x and then with respect to y. =
=A0(I would
> > not have used the term "cumulative probability" to describe this particu=
lar
> > probability, however.) =A0Unfortunately, matlab's 'dblquad' is unable to=

> evaluate
> > a double integral in this form where the inner limits of integration dep=
end
> on
> > the outer variable.
>
> > =A0 I can see two possibilities for numerically calculating this using m=
atlab. =A0
> > First, if you already know the cumulative distribution of X as a functio=
n, F
> (x),
> > the above can then be directly expressed as the single integral:
>
> > =A0integral, 0 to inf, (F(y+k)-F(y))*g(y) dy
>
> > by direct substitution of F(y+k)-F(y) for the inner integral after facto=
ring
> out
> > the g(y). =A0Any of the matlab single quadrature functions can be used f=
or
> this:
> > 'quad', 'quad8', etc.
>
> > =A0 The other possibility is to make a change of variable. =A0Define
>
> > =A0u =3D x - y
>
> > Then the iterated integral can be converted to the form
>
> > =A0integral, 0 to inf, integral, 0 to k, f(u+y)*g(y) du dy
>
> > (Note that the Jacobian in this case is just 1.) =A0This is in a form th=
at can be
> > used by 'dblquad' since the inner limits are independent of y values. =
=A0This is
> > the method you should use if you know only the two probability density
> > functions.
>
> > =A0 Note that each of these forms requires a separate numerical integrat=
ion
> for
> > each desired value of k. =A0There is always the possibility that, depend=
ing on
> > the functions, f(x) and g(y), there may exist an analytic solution to th=
ese
> > integrals using the 'int' function of the Symbolic Toolbox, but the odds=
 are
> > heavily against this I would think.
>
> > Roger Stafford
>
> -------
> =A0 I should have warned you that matlab's numerical quadrature functions
> cannot actually accept infinite limits. =A0in your problem you will have t=
o use a
> finite limit that is sufficiently large to encompass all the areas with si=
gnificant
> probability densities, but not so large as to confuse the routine into tak=
ing
> too few samples of the integrands in the critical areas to attain the need=
ed
> accuracy. =A0These routines don't seem to be particularly robust in this r=
egard. =A0
> You'll have to experiment with them a bit.
>
> Roger Stafford- Hide quoted text -
>
> - Show quoted text -

Another way is to hit it with the sledge hammer: Monte Carlo
simulation.
Sample, say, a million from each distribution, then do a 2-D histogram
on the results.
Repeat this many times and use allstats (from the File Exchange) to
calculate the statistics for each bin in the histogram.
When the standard errors for each bin in the histogram reduce to an
acceptable level, you're done.