[R] Poisson-lognormal probability calculations

[k.jewell at campden.co.uk]

Hi,

just for the record, although I don't think it's relevant (!)
-------------------------------------
> sessionInfo()
R version 2.6.0 (2007-10-03)
i386-pc-mingw32

locale:
LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United 
Kingdom.1252;LC_MONETARY=English_United 
Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252

attached base packages:
[1] stats4    splines   stats     graphics  grDevices utils     datasets 
methods   base

other attached packages:
[1] VGAM_0.7-5         xlsReadWrite_1.3.2
--------------------------------

I'm having some problems with a Poisson-lognormal density (mass?) function.

VGAM has the dpolono function, but that doesn't work for x-values over 170,
and I need to go to *much* bigger numbers. It fails first because of gamma
overflow, then because of non-finite integrand.
---------------------
> VGAM::dpolono(170)
[1] 4.808781e-09
> VGAM::dpolono(171)
[1] 0
Warning message:
In VGAM::dpolono(171) : value out of range in 'gammafn'
> VGAM::dpolono(172)
[1] 0
Warning message:
In VGAM::dpolono(172) : value out of range in 'gammafn'
> VGAM::dpolono(173)
[1] 0
Warning message:
In VGAM::dpolono(173) : value out of range in 'gammafn'
> VGAM::dpolono(174)
Error in integrate(f = integrand, lower = -Inf, upper = Inf, x = x[i],  :
  non-finite function value
-------------

I tidied up a little (to my eyes only, no offence intended to the estimable
authors of VGAM) and avoided the gamma overflow by using logs - OK so far,
it agrees almost perfectly with VGAM::dpolono for x up to 170, and now
extends the range up to 173 (wow!!).
-------------
dApolono <-
function (x, meanlog = 0, sdlog = 1, ...) {
    require(stats)
    integrand <- function(t, x, meanlog, sdlog)
exp(-t+x*log(t)-log(sdlog*t*sqrt(2*pi))-0.5*((log(t)-meanlog)/sdlog)^2)
    mapply(function(x, meanlog, sdlog, ...){
       temp <- try(
         integrate(f = integrand, lower = 0, upper = Inf, x = x, meanlog =
meanlog, sdlog = sdlog, ...)
                  )
         ifelse(inherits(temp, "try-error"), NA,
exp(log(temp$value)-lgamma(x+1)))
                                      }
              , x, meanlog, sdlog, ...
           )
                                           }
plot(log(dApolono(0:173)))
dApolono(174)
-----------------------------

Addressing the non-finite integrand, I noticed that the gamma(x+1) divisor
was outside the integrand so that as x gets bigger the integrand gets
bigger, only to be finally divided by an increasing gamma(x+1). I reasoned
that putting the divisor inside the integral might incur a substantial
performance hit, but would keep the integrand at reasonable values.

--------------------------------------
dApolono <-
function (x, meanlog = 0, sdlog = 1, ...) {
    require(stats)
    integrand <- function(t, x, meanlog, sdlog)
exp(-t+x*log(t)-log(sdlog*t*sqrt(2*pi))-0.5*((log(t)-meanlog)/sdlog)^2-lgamm
a(x+1))
    mapply(function(x, meanlog, sdlog, ...){
       temp <- try(
         integrate(f = integrand, lower = 0, upper = Inf, x = x, meanlog =
meanlog, sdlog = sdlog, ...)
                  )
         ifelse(inherits(temp, "try-error"), NA, temp$value)
                                      }
              , x, meanlog, sdlog, ...
           )
                                           }
plot(log(dApolono(0:173)))
dApolono(174)
dApolono(1E3)
-------------------------

This avoids the non-finite integrand and gives me answers at much higher x,
but now at least some of the answers are wrong (not to say silly), even in
the range where the other versions worked. I've tried other variations
including changing the variable of integration to log(t) and integrating
dpois()*dlnorm(), but I can't fix it; if the factorial (=gamma) is inside
the integrand I get silly answers, if it's outside I get non-finite
integrand.

I'm tearing my hair. Can anyone suggest where I may be going wrong? Any
suggestions at all will be appreciated.

Thanks in advance,

Keith Jewell
mailto:k.jewell at campden.co.uk telephone (direct) +44 (0)1386 842055




_____________________________________________________________________
The information in this document and attachments is given after the exercise of all reasonable care and skill in its compilation, preparation and issue, but is provided without liability in its application or use.  It may contain privileged information that is exempt from disclosure by law and may be confidential.  If you are not the intended recipient you must not copy, distribute or take any action in reliance on it.  If you have received this document in error please notify us and delete this message from your system immediately.

Campden & Chorleywood Food Research Association Group;
Campden & Chorleywood Food Research Association (company limited by guarantee),registered number 510618);
CCFRA Group Services Ltd. (registered number 3841905); and
CCFRA Technology Ltd  (registered number 3836922),
all registered in England and Wales with the registered office at Station Road, Chipping Campden, Gloucestershire, GL55 6LD.

The Group may monitor e-mail traffic data and also the content of e-mail for the purposes of security and staff training.

Any advice given is subject to our normal terms and conditions of trading, a copy of which is available on request.

________________________________________________________________________
This e-mail has been scanned for all viruses by MessageL...{{dropped:2}}