# [R] Fitting a Bivariate Poisson log-normal distribution

Famoye, Felix famoy1kf at cmich.edu
Sun Jul 27 21:36:41 CEST 2008

```I wrote an R program to fit a bivariate Poisson log-normal model. The
model is first proposed by Aitchison and Ho (1989) in Biometrika, Volume
76 pages 643-653. The model involves using a double integration. The
following is the program and the data that I used. My major problem was
that R was running for a long time (more than 3 hours) with no results.
I would like to place some printing at some places to see if R was doing
something. Unfortunately, there was nothing printed when I included
printing in the function "f". Note that I have used the package "adapt"
for double integration. Can someone provide some help on this problem? I
would like to know if R is doing anything at all in fitting the model.
Thank you for your time.

#This program will be used to estimate the parameters of
Poisson-lognormal model.

y1 = yd[[1]]
y2 = yd[[2]]
n = length(y1)
x0= rep(1, n)
y1  = cbind(y1)
y2  = cbind(y2)
yy=cbind(y1,y2)
xx = cbind(x0,yd[[3]])
av1=mean(y1)
av2=mean(y2)
cov=var(yy)
mu=cbind(av1, av2)
rho=cov[1,2]/sqrt(cov[1,1]*cov[2,2])
sd1=sqrt(cov[1,1])
sd2=sqrt(cov[2,2])
oth=rbind(sd1,sd2,rho)
z = vector(length=2)  #This is a column vector with two rows

#To find the inital estimates for the regression coefficients
para1 = ginv(t(xx) %*% xx) %*% (t(xx) %*% log(y1+0.5))
para2 = ginv(t(xx) %*% xx) %*% (t(xx) %*% log(y2+0.5))
parab = rbind(para1, para2, oth)
parab

#---------------------------------------------------------------------
# Fitting Poisson lognorma model
#f is the function that defines the (negative) log likelihood
f <- function(parab,y1,y2,xx)
{
vv1=(parab[5])^2
vv2=(parab[6])^2
para.1= cbind(parab[1:2])
para.2= cbind(parab[3:4])
mu.1 = xx%*%para.1 - vv1/2
#mu.1 = exp(mu.1)
mu.2 = xx%*%para.2 - vv2/2
#mu.2 = exp(mu.2)
mu.v=cbind(mu.1,mu.2)
cov[1,1]=vv1
cov[2,2]=vv2
cov[1,2]=parab[7]*parab[5]*parab[6]
cov[2,1]=cov[1,2]
n=length(y1)
va=rep(0.0, n)
for (i in 1:n)
{
fred=function(z)
{
tz=log(z)-mu.v[i,]
qq=(t(tz)%*%ginv(cov))%*%tz
fq=(exp(-qq/2))/(2*pi*z[1]*z[2]*sqrt(det(cov)))
p1=((z[1]^y1[i])*exp(-z[1]))/gamma(y1[i]+1)
p2=((z[2]^y2[i])*exp(-z[2]))/gamma(y2[i]+1)
int=fq*p1*p2
}
v=mu.v,cov=cov,y1=y1,y2=y2)\$value
}
-sum(va)
}

I.p = nlminb(start=parab,objective=f,y1=y1,y2=y2,xx=xx)
print("Poisson Log-normal Regression Model")
I.p

0   0   0
0   0   0
0   0   0
0   0   0
0   0   0
0   0   0
0   0   0
0   0   0
0   0   0
0   0   0
0   0   0
0   0   0
1   0   0
1   0   0
1   0   0
1   0   0
1   0   0
1   0   0
1   0   0
2   0   0
2   0   0
2   0   0
2   0   0
2   0   0
2   0   0
3   0   0
3   0   0
3   0   0
3   0   0
3   0   0
3   0   0
4   0   0
4   0   0
4   0   0
5   0   0
5   0   0
7   0   0
8   0   0
9   0   0
9   0   0
12   0   0
13   0   0
14   0   0
16   0   0
20   0   0
0   1   0
0   1   0
0   1   0
0   1   0
0   1   0
0   1   0
0   1   0
1   1   0
1   1   0
1   1   0
1   1   0
1   1   0
1   1   0
1   1   0
1   1   0
1   1   0
1   1   0
1   1   0
1   1   0
1   1   0
2   1   0
2   1   0
2   1   0
2   1   0
2   1   0
2   1   0
2   1   0
3   1   0
3   1   0
3   1   0
4   1   0
4   1   0
5   1   0
5   1   0
5   1   0
9   1   0
10   1   0
10   1   0
15   1   0
0   2   0
0   2   0
0   2   0
0   2   0
0   2   0
0   2   0
0   2   0
0   2   0
0   2   0
1   2   0
1   2   0
1   2   0
1   2   0
1   2   0
1   2   0
1   2   0
1   2   0
1   2   0
2   2   0
2   2   0
2   2   0
3   2   0
4   2   0
4   2   0
5   2   0
6   2   0
6   2   0
6   2   0
7   2   0
7   2   0
12   2   0
15   2   0
26   2   0
0   3   0
0   3   0
0   3   0
0   3   0
0   3   0
0   3   0
0   3   0
0   3   0
0   3   0
1   3   0
1   3   0
1   3   0
1   3   0
1   3   0
2   3   0
2   3   0
4   3   0
5   3   0
7   3   0
10   3   0
0   4   0
0   4   0
0   4   0
0   4   0
1   4   0
1   4   0
1   4   0
1   4   0
3   4   0
3   4   0
0   5   0
0   6   0
0   6   0
1   6   0
3   6   0
0   7   0
0   7   0
1   7   0
0   8   0
1   8   0
0   9   0
0   9   0
0   9   0
0   0   1
0   0   1
0   0   1
0   0   1
0   0   1
0   0   1
0   0   1
0   0   1
1   0   1
1   0   1
1   0   1
1   0   1
1   0   1
1   0   1
1   0   1
1   0   1
1   0   1
1   0   1
1   0   1
1   0   1
2   0   1
2   0   1
2   0   1
2   0   1
2   0   1
2   0   1
3   0   1
3   0   1
3   0   1
3   0   1
3   0   1
3   0   1
3   0   1
3   0   1
4   0   1
4   0   1
4   0   1
4   0   1
4   0   1
4   0   1
4   0   1
4   0   1
5   0   1
5   0   1
7   0   1
7   0   1
8   0   1
9   0   1
9   0   1
11   0   1
15   0   1
15   0   1
19   0   1
20   0   1
0   1   1
0   1   1
0   1   1
0   1   1
0   1   1
0   1   1
0   1   1
0   1   1
0   1   1
0   1   1
0   1   1
1   1   1
1   1   1
1   1   1
2   1   1
2   1   1
2   1   1
2   1   1
2   1   1
4   1   1
4   1   1
5   1   1
5   1   1
5   1   1
6   1   1
8   1   1
8   1   1
8   1   1
8   1   1
9   1   1
23   1   1
23   1   1
0   2   1
0   2   1
0   2   1
0   2   1
0   2   1
0   2   1
0   2   1
0   2   1
0   2   1
0   2   1
0   2   1
0   2   1
0   2   1
0   2   1
0   2   1
1   2   1
1   2   1
1   2   1
1   2   1
1   2   1
1   2   1
1   2   1
1   2   1
1   2   1
1   2   1
2   2   1
2   2   1
2   2   1
3   2   1
3   2   1
3   2   1
4   2   1
4   2   1
6   2   1
0   3   1
0   3   1
0   3   1
0   3   1
0   3   1
0   3   1
0   3   1
1   3   1
1   3   1
1   3   1
1   3   1
1   3   1
2   3   1
2   3   1
2   3   1
3   3   1
6   3   1
0   4   1
0   4   1
0   4   1
0   4   1
0   4   1
1   4   1
1   4   1
1   4   1
1   4   1
2   4   1
5   4   1
0   5   1
0   5   1
0   5   1
1   5   1
1   5   1
0   6   1
0   6   1
1   7   1
1   7   1
1   8   1
0  10   1
0  10   1

-----
------------------------------------------------------
Felix Famoye
Department of Mathematics
Central Michigan University
Mt. Pleasant, Michigan 48859-0001

e-mail:    felix.famoye at cmich.edu
web site:  http://www.cst.cmich.edu/users/famoy1kf/
voice:     (989)774-5497,   fax: (989)774-2414

```