[R] Erlang Loss Function in R

[Tom Potter]

This is my first post here.  As a new user of R, I was pleased to discover a
simple and  elegant way to compute in R the Erlang Loss Function, which I
present below.

INTRODUCTION
The Erlang Loss Function is defined as

elf(n, a ) = ( a^n/n! ) / ( 1 + a + a^2/2 + a^3/3! + ... + a^n/n!)

for all real numbers a>=0 and all integers n>=0 (except that the function is
undefined when  BOTH n and a are zero).  The function can be extended
analytically to non-integral values  of n by various means.  The basic
reference is: Jagerman, D. L. (1974), "Some properties of  the Erlang loss
function", Bell System Tech. J. (53), 525-551.

The Erlang Loss Function is important in queueing theory, especially in
telecom  applications.  One use is in computing the blockage in the M/M/c/c
queue (Poisson arrivals,  exponential service times, finite number of
servers, no waiting positions).  Here

  blockage = elf( nserver, load )

where

  nserver = number of servers
  load = offered load in erlangs
  blockage = probability that an arriving customer finds all servers busy

FIRST ATTEMPT: ITERATION
Using a familiar iteration ("Erlang's B formula") one could write this code:

elf <- function(n, a)
  {
    b <- 1 ;  i<- 0
    while((i <- i+1)<=n) b <- b*a/(i+b*a)
    b
  }

This works for all integers n>=0 and all real numbers a>=0.  But it has
these deficiencies:  (1) slow for large 'n', (2) gives incorrect values when
'n' is not an integer, (3) will not  work when 'n' is a vector (but it does
allow 'a' to be a vector).

CODING USING GAMMA DISTRIBUTION
A nicer version of elf() exploits R's wonderful gamma distribution
functions!  I use the  fact, perhaps not generally known, that the Erlang
Loss Function is the ratio of the gamma  density to the gamma CDF, for
suitable parameters.  (Also, I'll change the name from elf()  to erlb(),
because this function is sometimes called "Erlang's B function", and because
'elf' is an abbreviation in frequent use with a different meaning in the
Unix world.)  Here  is the code:

erlb<- function(nsrv, nerl, log=FALSE)
{
  # Returns Erlang-B blockage for nsrv servers and nerl erlangs of traffic.
  # If log=TRUE, then logarithm of the blockage is returned.  Works for all
  # nsrv>=0 and all nerl>=0.  Both nsrv and nerl can be non-integers.

  lnerlb<-
    dgamma(nerl, shape=nsrv+1, scale=1, log=TRUE) -
    pgamma(nerl, shape=nsrv+1, scale=1, log.p=TRUE, lower.tail=FALSE) ;
  if (log) lnerlb else exp(lnerlb)
}

This version has these advantages
.  MUCH faster, especially for large values of the 1st argument
.  Works for non-integral 'nsrv'
.  Works when either 'nsrv', or 'nerl', or both are vectors
.  Can return an accurate logarithm of the Erlang Loss Function if that is
needed.

Tom Potter
Nashville, Tennessee
tpotter at home.com <mailto:tpotter at home.com>

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._