Prev: clustering 'locally' for large data set?
Next: restore the image curvation
From: Bradley on 25 May 2010 19:05
Hi guys,

I'm trying to compute the below double integral with variable limits using quad2d;


quad2d(@(u,v) (1./u).*lognpdf(1./u,mu_a,sigma_a).*(1./(1-logncdf(((1/mc).*(delta.*(1./v)+c*(eta/(1-eta))*chi*theta-(c/((1-eta)*q)))),mu_a,sigma_a))).*lognpdf(1./v,mu_h,sigma_h).*(-(1./u).^2).*(-(1./v).^2).*(1/(1-logncdf(hbar,mu_h,sigma_h))),(1./hbar),0,@(v)1./((1/mc).*(delta.*(1./v)+c*(eta/(1-eta))*chi*theta-(c/((1-eta)*q)))),0)

where u and v are the variables that are integrated over and everything else are constants, given by;

mu_a = 0;
sigma_a = 0.1;
mc = 0.9091;
delta = 0.8510;
c = 0.049;
eta = 0.5;
chi = 0.17;
theta = 0.8571;
q = 0.7;
mu_h = 0;
sigma_h = 0.1;
hbar = 1.0707;

quad2d will not be able to compute the above and reports that a singularity is likely.
Does anybody know how to work around such singularities?

Many thanks,
Brad.
From: Bradley on 25 May 2010 19:34
"Bradley " <bspeigner_1(a)hotmail.com> wrote in message <hthl34$gsn$1(a)fred.mathworks.com>...
> Hi guys,
>
> I'm trying to compute the below double integral with variable limits using quad2d;
>
>
> quad2d(@(u,v) (1./u).*lognpdf(1./u,mu_a,sigma_a).*(1./(1-logncdf(((1/mc).*(delta.*(1./v)+c*(eta/(1-eta))*chi*theta-(c/((1-eta)*q)))),mu_a,sigma_a))).*lognpdf(1./v,mu_h,sigma_h).*(-(1./u).^2).*(-(1./v).^2).*(1/(1-logncdf(hbar,mu_h,sigma_h))),(1./hbar),0,@(v)1./((1/mc).*(delta.*(1./v)+c*(eta/(1-eta))*chi*theta-(c/((1-eta)*q)))),0)
>
> where u and v are the variables that are integrated over and everything else are constants, given by;
>
> mu_a = 0;
> sigma_a = 0.1;
> mc = 0.9091;
> delta = 0.8510;
> c = 0.049;
> eta = 0.5;
> chi = 0.17;
> theta = 0.8571;
> q = 0.7;
> mu_h = 0;
> sigma_h = 0.1;
> hbar = 1.0707;
>
> quad2d will not be able to compute the above and reports that a singularity is likely.
> Does anybody know how to work around such singularities?
>
> Many thanks,
> Brad.

Forgot to note that I'm looking for workarounds that do not involve changing the limits. So, for example, changing the zero limits in the above integral to a small positive number would get rid of the problem.

It's been a long day/month/year/etc...
From: Walter Roberson on 25 May 2010 22:33
Bradley wrote:

> I'm trying to compute the below double integral with variable limits
> using quad2d;

> quad2d(@(u,v)
> (1./u).*lognpdf(1./u,mu_a,sigma_a).*(1./(1-logncdf(((1/mc).*(delta.*(1./v)+c*(eta/(1-eta))*chi*theta-(c/((1-eta)*q)))),mu_a,sigma_a))).*lognpdf(1./v,mu_h,sigma_h).*(-(1./u).^2).*(-(1./v).^2).*(1/(1-logncdf(hbar,mu_h,sigma_h))),(1./hbar),0,@(v)1./((1/mc).*(delta.*(1./v)+c*(eta/(1-eta))*chi*theta-(c/((1-eta)*q)))),0)
>
>
> where u and v are the variables that are integrated over and everything
> else are constants

> quad2d will not be able to compute the above and reports that a
> singularity is likely.

If I have parsed the documentation properly,
http://www.mathworks.co.uk/access/helpdesk/help/toolbox/stats/lognpdf.html

then lognpdf(1/u, mu_a, sigma_a) should be the same as
lognpdf(u, -mu_a, sigma_a)

The only term involve x is also the only term involving mu, and is
(ln x - mu)^2 . When you put in x = 1/u then ln 1/u = -ln u, so that
would become (-ln u - mu)^2 which would be ((-1)*(ln u + mu))^2 which is
(ln u + mu)^2 -- which is the same thing as if you had put in u for x
originally and had put in the negative of mu.

You can thus eliminate the singularities for 1/u and 1/v in the lognpdf
calls.

However, the leading (1/u) term would combine with the (-(1./u).^2) term
of the denominator to give a 1/u^3 term, and the (-(1./v).^2) term of
the denominator is likewise not canceled out by anything. In order for
you to be able to cancel the singularity due to the 1/u^3, you would
have to show that lognpdf(1/u, mu_a, sigma_a) has u in the numerator to
a higher power than u^3. Unfortunately it does not: instead the
expression would be a constant divided by u^5. The singularity due to u
becoming 0 is thus unavoidable, and likewise the singularity due to v
becoming 0 is unavoidable (it would involve constant/v^4)
 | 
Pages: 1
Prev: clustering 'locally' for large data set?
Next: restore the image curvation