if (FALSE) {
## Means of a BGWM process 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 )
# mean matrix of the process
I.matriz.m <- BGWM.mean(Dists.i, "independents", d)
# mean vector of the population in the nth generation
# from vector N representing the initial population
I.vector.m.n_N <- BGWM.mean(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 )
# mean matrix of the process
M.matrix.m <- BGWM.mean(Dists.m, "multinomial", d)
# mean vector of the population in the nth generation
# from vector N representing the initial population
M.vector.m.n_N <- BGWM.mean(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 )
# mean matrix of the process
G.matrix.m <- BGWM.mean(Dists.g, "general", d)
# mean vector of the population in the nth generation
# from vector N representing the initial population
G.vector.m.n_N <- BGWM.mean(Dists.g, "general", d, n, N)
# Comparison of results
I.vector.m.n_N
I.vector.m.n_N - M.vector.m.n_N
M.vector.m.n_N - G.vector.m.n_N
G.vector.m.n_N - I.vector.m.n_N
}
Run the code above in your browser using DataLab