rBGWM {Branching}R Documentation

Simulating a multi-type Bienayme - Galton - Watson process

Description

Generate the trajectories of a multi-type Bienayme - Galton - Watson process from its offspring distributions, using three different algorithms based on three different classes or families of these processes.

Usage

rBGWM(dists, type=c("general","multinomial","independents"), d,
      n, z0=rep(1,d), c.s=TRUE, tt.s=TRUE, rf.s=TRUE, file=NULL)

Arguments

dists

offspring distributions. Its structure depends on the class of the Bienayme - Galton - Watson process (See details and examples).

type

Class or family of the Bienayme - Galton - Watson process (See details).

d

positive integer, number of types.

n

positive integer, maximum lenght of the wanted trajectory.

z0

nonnegative integer vector of size d; initial population by type.

c.s

logical value, if TRUE, the output object will include the generated trajectory of the process with the number of individuals for every combination parent type - descendent type.

tt.s

logical value, if TRUE, the output object will include the generated trajectory of the process with the number of descendents by type.

rf.s

logical value, if TRUE, the output object will include the generated trajectory of the process with the relative frequencies by type.

file

the name of the output file where the generated trajectory of the process with the number of individuals for every combination parent type - descendent type could be stored.

Details

This function performs a simulation of a multi-type Bienayme - Galton - Watson process (BGWM) from its offspring distributions.

From particular offspring distributions and taking into account a differentiated algorithmic approach, we propose the following classes or types for these processes:

general This option is for BGWM processes without conditions over the offspring distributions, in this case, it is required as input data for each distribution, all d-dimensional vectors with their respective, greater than zero, probability.

multinomial This option is for BGMW processes where each offspring distribution is a multinomial distribution with a random number of trials, in this case, it is required as input data, d univariate distributions related to the random number of trials for each multinomial distribution and a d \times d matrix where each row contains probabilities of the d possible outcomes for each multinomial distribution.

independents This option is for BGMW processes where each offspring distribution is a joint distribution of d combined independent discrete random variables, one for each type of individuals, in this case, it is required as input data d^2 univariate distributions.

The structure need it for each classification is illustrated in the examples.

These are the univariate distributions available:

unif Discrete uniform distribution, parameters min and max. All the non-negative integers between min y max have the same probability.

binom Binomial distribution, parameters n and p.

p(x) = choose(n,x) p^x (1-p)^(n-x)

for x = 0, ..., n.

hyper Hypergeometric distribution, parameters m (the number of white balls in the urn), n (the number of white balls in the urn), k (the number of balls drawn from the urn).

p(x) = choose(m, x) choose(n, k-x) / choose(m+n, k)

for x = 0, ..., k.

geom Geometric distribution, parameter p.

p(x) = p (1-p)^x

for x = 0, 1, 2, ...

nbinom Negative binomial distribution, parameters n and p.

p(x) = Gamma(x+n)/(Gamma(n) x!) p^n (1-p)^x

for x = 0, 1, 2, ...

pois Poisson distribution, parameter lambda.

p(x) = lambda^x exp(-lambda)/x!

for x = 0, 1, 2, ...

norm Normal distribution rounded to integer values and negative values become 0, parameters mu and sigma.

p(x) = \int_{x-0.5}^{x+0.5} 1/(sqrt(2 pi) sigma) e^-((t - mu)^2/(2sigma^2)) dt

for x = 1, 2, ...

p(x) = \int_{-∞}^{0.5} 1/(sqrt(2 pi) sigma) e^-((t - mu)^2/(2sigma^2)) dt

for x = 0

lnorm Lognormal distribution rounded to integer values, parameters logmean = mu y logsd = sigma.

p(x) = \int_{x-0.5}^{x+0.5} 1/(sqrt(2 pi) sigma t) e^-((log t - mu)^2 / (2sigma^2))dt

for x = 1, 2, ...

p(x) = \int_{0}^{0.5} 1/(sqrt(2 pi) sigma t) e^-((log t - mu)^2 / (2sigma^2)) dt

for x = 0

gamma Gamma distribution rounded to integer values, parameters shape = a y scale = s.

p(x)= \int_{x-0.5}^{x+0.5} 1/(s^a Gamma(a)) t^(a-1) e^-(t/s) dt

para x = 1, 2, ...

p(x)= \int_{0}^{0.5} 1/(s^a Gamma(a)) t^(a-1) e^-(t/s) dt

for x = 0

Value

An object of class list with these components:

i.d

input. number of types.

i.dists

input. offspring distributions.

i.n

input. maximum lenght of the generated trajectory.

i.z0

input. initial population by type.

o.c.s

output. A matrix indicating the generated trajectory of the process with the number of individuals for every combination parent type - descendent type.

o.tt.s

output. A matrix indicating the generated trajectory of the process with the number of descendents by type.

o.rf.s

output. A matrix indicating the generated trajectory of the process with the relative frequencies by type.

Author(s)

Camilo Jose Torres-Jimenez cjtorresj@unal.edu.co

References

Torres-Jimenez, C. J. (2010), Relative frequencies and parameter estimation in multi-type Bienaym? - Galton - Watson processes, Master's Thesis, Master of Science in Statistics. Universidad Nacional de Colombia. Bogota, Colombia.

Stefanescu, C. (1998), 'Simulation of a multitype Galton-Watson chain', Simulation Practice and Theory 6(7), 657-663.

Athreya, K. & Ney, P. (1972), Branching Processes, Springer-Verlag.

See Also

BGWM.mean, BGWM.covar, BGWM.mean.estim, BGWM.covar.estim

Examples

## Not run: 
## Simulation based on a model analyzed in Stefanescu(1998)

# Variables and parameters
d <- 2
n <- 30
N <- c(90, 10)
a <- c(0.2, 0.3)

# with independent distributions
Dists.i <- data.frame( name=rep( "pois", d*d ),
                       param1=rep( a, rep(d,d) ),
                       stringsAsFactors=FALSE )

rA <- rBGWM(Dists.i, "independents", d, n, N)

# with multinomial distributions
dist <- data.frame( name=rep( "pois", d ),
                    param1=a*d,
                    stringsAsFactors=FALSE )
matrix.b <- matrix( rep(0.5, 4), nrow=2 )
Dists.m <- list( dists.eta=dist, matrix.B=matrix.b )

rB <- rBGWM(Dists.m, "multinomial", d, n, N)

# with general distributions (approximation)
max <- 30
A <- t(expand.grid(c(0:max),c(0:max)))
aux1 <- factorial(A)
aux1 <- apply(aux1,2,prod)
aux2 <- apply(A,2,sum)
distp <- function(x,y,z){ exp(-d*x)*(x^y)/z }
p <- sapply( a, distp, aux2, aux1 )
prob <- list( dist1=p[,1], dist2=p[,2] )
size <- list( dist1=ncol(A), dist2=ncol(A) )
vect <- list( dist1=t(A), dist2=t(A) )
Dists.g <- list( sizes=size, probs=prob, vects=vect )

rC <- rBGWM(Dists.g, "general", d, n, N)

# Comparison chart
dev.new()
plot.ts(rA$o.tt.s,main="with independents")
dev.new()
plot.ts(rB$o.tt.s,main="with multinomial")
dev.new()
plot.ts(rC$o.tt.s,main="with general (aprox.)")

## End(Not run)

[Package Branching version 0.9.4 Index]