[R] a prolem with constrOptim

[Ravi Varadhan]

Hi,

I noticed that you are setting very small tolerance as convergence
criterion. In most optimization problems, it does not make much practical
sense to set such a low tolerance.  If you increase this, then `constrOptim'
also works. For example,

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

However, it should be noted that my function `contrOptim.nl' provides a
better answer for both data1 and data2 in terms of having a larger
log-likelihood at convergence.

Ravi.

----------------------------------------------------------------------------
-------

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvaradhan at jhmi.edu

Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
tml

 

----------------------------------------------------------------------------
--------


-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of Davidov, Ori (NIH/NIEHS) [V]
Sent: Monday, November 02, 2009 11:36 AM
To: r-help at r-project.org
Subject: [R] a prolem with constrOptim

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 at 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.