Famoye, Felix
2008-Jul-27 19:36 UTC
[R] Fitting a Bivariate Poisson log-normal distribution
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.
yd = read.table(file="H:\\Bivariate\\bact.txt")
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
}
va[i]=adapt(2,lo=c(0,0),up=c(100,100),functn=fred,min=1000,eps=1.e-6,mu.
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