| 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) = {n \choose x} {p}^{x} {(1-p)}^{n-x}
for x = 0, \dots, 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) = \left. {m \choose x}{n \choose k-x} \right/ {m+n \choose k}%
for x = 0, ..., k.
geom Geometric distribution, parameter p.
p(x) = p {(1-p)}^{x}
for x = 0, 1, 2, \dots
nbinom Negative binomial distribution, parameters n and
p.
p(x) = \frac{\Gamma(x+n)}{\Gamma(n) x!} p^n (1-p)^x
for x = 0, 1, 2, \dots
pois Poisson distribution, parameter \lambda.
p(x) = \frac{\lambda^x e^{-\lambda}}{x!}
for x = 0, 1, 2, \dots
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} \frac{1}{\sqrt{2\pi}\sigma} e^{-(t-\mu)^2/2\sigma^2}dt%
for x = 1, 2, \dots
p(x) = \int_{-\infty}^{0.5} \frac{1}{\sqrt{2\pi}\sigma} e^{-(t-\mu)^2/2\sigma^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} \frac{1}{\sqrt{2\pi}\sigma t} e^{-(\log(t) - \mu)^2/2 \sigma^2 }dt%
for x = 1, 2, \dots
p(x) = \int_{0}^{0.5} \frac{1}{\sqrt{2\pi}\sigma t} e^{-(\log(t) - \mu)^2/2 \sigma^2 }dt%
for x = 0
gamma Gamma distribution rounded to integer values,
parameters shape =\alpha y scale =\sigma.
p(x)= \int_{x-0.5}^{x+0.5} \frac{1}{{\sigma}^{\alpha}\Gamma(\alpha)}{t}^{\alpha-1} e^{-t/\sigma} dt%
para x = 1, 2, \dots
p(x)= \int_{0}^{0.5} \frac{1}{{\sigma}^{\alpha}\Gamma(\alpha)}{t}^{\alpha-1} e^{-t/\sigma} 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 Bienayme - 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)