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: rvarad...@jhmi.edu
Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
tml
----------------------------------------------------------------------------
--------
-----Original Message-----
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Davidov, Ori (NIH/NIEHS) [V]
Sent: Monday, November 02, 2009 11:36 AM
To: r-help@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@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.
______________________________________________
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.
