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