Hi,

I apologize for the long message but the problem I encountered can't be stated in a few lines.

I am having some problems with the function constrOptim. My goal is to maximize the likelihood of product of K multinomials, each with four catagories under linear constraints on the parameter values. I have found that the function does not work for many data configurations.

#The likelihood and score functions are:

likelihood = function(theta,data)
{
                p = matrix(theta,nrow=3,byrow=F)
                s = 1-apply(p,2,sum)
                p = rbind(p,s)
                Z = matrix(data,nrow=4,byrow=F)
                sum(log(p)*Z)
}

score = function(theta,data)
{
                k = length(data)/4
                S = numeric(3*k)
                for(i in 1:k)
                {
                                n = data[(1+4*(i-1)):(4*i)]
                                p = theta[(1+3*(i-1)):(3*i)]
                                P = 1-sum(p)
                                S[(1+3*(i-1)):(3*i)] = n[1:3]/p-n[4]/P
                }
                S
}

#where theta=(p11,p12,p13,p21,p22,p23,...,pK1,pK2,pK3).

#The function Rmat calculates the restriction matrix needed for constrained estimation

Rmat = function(k)
{
                R = matrix(1,4,3)
                R[1,2] = R[1,3] = R[2,3] = R[3,2] = 0
                RR = cbind(-R,R)
                RRR = matrix(0,4*(k-1),3*k)
                for ( i in 1:(k-1) ) RRR[4*(i-1)+1:4,3*(i-1)+1:6] = RR
                RRR
}

#The function gen.data generates data

gen.data = function(p,N)
{
                k = ncol(p)
                Data = matrix(0,4,k)
                for(j in 1:k)
                {
                                Data[,j] = as.vector(rmultinom(1, size = N[j], prob=p[,j]))
                }
                as.vector(Data)
}

#I have found that the function returns an error under some configurations of the data ( 0's seem to harm it). Why is this happening and how can it #be corrected? Why do some configurations with 0's crash the program and others don't?

# For example for the data:

data1 = c(0,6,8,6,3,4,4,9,2,2,5,11)
data2 = c(0,6,8,6,0,4,4,9,2,2,5,11)

# inputs for constrOptim

K = 3
RR = Rmat(K)
zeros = rep(0,nrow(RR))
theta.0 = numeric(3*K)
for(i in 1:K) {theta.0[(1+3*(i-1)):(3*i)] = rep(1/10,3)+i/100}

# data1 gives an error but data2 does not!

constrOptim(theta=theta.0, f=likelihood, grad=score, ui=RR,ci=zeros,
                data = data1, method="BFGS",
                control=list(fnscale=-1,abstol=1e-16,reltol=1e-16))$par

constrOptim(theta=theta.0, f=likelihood, grad=score, ui=RR,ci=zeros,
                data = data2, method="BFGS",
                control=list(fnscale=-1,abstol=1e-16,reltol=1e-16))$par


# If you want to further check run the simulation below.

K = 3
N = rep(20,K)
p = c(1/10,2/10,2/10,5/10)
p = matrix(rep(p,K),nrow=4,byrow=F)


RR = Rmat(K)
zeros = rep(0,nrow(RR))
theta.0 = numeric(3*K)
for(i in 1:K) {theta.0[(1+3*(i-1)):(3*i)] = rep(1/10,3)+i/100}

for(i in 1:1000)
{
                data = gen.data(p,N)
                print(i)
                print(data)
                theta.til = constrOptim(theta=theta.0, f=likelihood, grad=score, ui=RR,ci=zeros,
                data = data, method="BFGS",
                control=list(fnscale=-1,abstol=1e-16,reltol=1e-16))$par

                print(theta.til)
}



	[[alternative HTML version deleted]]


