[R] a prolem with constrOptim

[Ravi Varadhan]

Hi Davidov,

You can use the `constrOptim.nl' function that I have written for solving
problems like yours.  It is better and more general than the `constrOptim'
function.  It is able to solve your problem, but I am not sure how well
posed your problem is.  One of the parameters seems to go to 0.  May be this
ok.

Here is the code for solving your problem:

	source("constrOptim_nl.txt")  
	
	hin <- function(par, ...) {
	# function that defines the inequality constraints
		   k = 3
                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
		   c(RRR %*% par)
	}

	ans1 <- constrOptim.nl(par=theta.0, fn=likelihood, hin=hin,
data=data1, control=list(fnscale=-1))

	ans2 <- constrOptim.nl(par=theta.0, fn=likelihood, hin=hin,
data=data2, control=list(fnscale=-1))

This works for both data1 and data2.

> ans1
$par
[1] 7.158690e-10 2.500022e-01 2.682024e-01 1.314431e-01 1.750843e-01
2.118290e-01 2.089612e-01 9.771505e-02 2.647332e-01

$value
[1] -72.81515


> ans2
$par
[1] 8.068277e-12 1.659718e-01 2.332808e-01 1.363536e-04 1.918888e-01
2.331728e-01 7.855021e-02 1.240325e-01 2.235095e-01

$value
[1] -65.75015

Let me know if you are interested in obtaining the `constrOptim.nl'
function.   

Hope this helps,
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.