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]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
