CODEWITHSUNDEEP
roptimization
dgranito10
dgranito10

Reputation: 3

Wrong result from constroptim function

I'm trying to use constrOptim to optimize the sum of square errors from a linear multiple regression. The main equation should be D = Beta1*Xa+Beta2*Xb+Beta3*Xc+Beta4*Xd , with D,Xa,Xb,Xc,Xd from a imported .csv file, and the Betas are the coefficients I want to find, minimizing the quadratic errors.

So far I imported the file.csv to R, named each column as Ds,Xa,Xb,Xc,Xd, created the objfunction= function(Beta1,Beta2,Beta3,Beta4)'sum(E²)'=(sum(D) - sum(Beta1*Xa+Beta2*Xb+Beta3*Xc+Beta4*Xd))^2)

created the matrix 'C' and vector 'd' to configure the constraints that should restrict the Beta's to <=0. I dont know how to find the feasible region, although I've used initial values that made the function work.

Here is the code:

> Tabela= read.table("Simulacao.csv", header=T, sep= ";")
> Tabela
   D A  B  C D.1
1 -1 1 -1  0   0
2  4 0  0  1  -1
3  4 1  0 -1   0
4  0 0  1  0  -1
5 -2 1  0  0  -1
> Ds= Tabela[,1]
> Xa= Tabela[,2]
> Xb= Tabela[,3]
> Xc= Tabela[,4]
> Xd= Tabela[,5]
> simulaf= function(x1,x2,x3,x4) {
+ Ds= Tabela[,1]
+ Xa= Tabela[,2]
+ Xb= Tabela[,3]
+ Xc= Tabela[,4]
+ Xd= Tabela[,5]
+ J=sum(Ds)
+ H=sum(x1*Xa+x2*Xb+x3*Xc+x4*Xd)
+ sx=(J-H)^2
+ return(sx)
+  }
> s= function(x)  {simulaf(x[1],x[2],x[3],x[4])}
> d= c(0,0,0,0)
> C= matrix(c(-1,0,0,0,0,-1,0,0,0,0,-1,0,0,0,0,-1),nrow=4,ncol=4,byrow=T)
> constrOptim(c(-1,-1,-1,-1),s,NULL,C,d)
$par
[1] -0.2608199 -0.8981110 -1.1095961 -1.9274866

The result I expect should be:

$par
[1] -0.125 0 -0.5 -0.875

After researching this, my conclusions are that it could be because I'm using bad initial values, parameterization problem (don't understand why its needed) or if it's simply that I have programmed it incorrectly.

What do I need to do to fix this?

Upvotes: 0

Views: 128

Answers (1)

Hong Ooi
Hong Ooi

Reputation: 57686

The formula for the sum of squared errors is

sum((y - yhat)^2)

and not

(sum(y) - sum(yhat))^2

where yhat is the predicted value.

Also, if your only constraints are that the estimated betas should be negative (which is a bit weird, usually you want them to be positive but never mind), then you don't need constrOptim. Regular optim(method="L-BFGS-B") or nlminb will work with so-called box constraints.

Upvotes: 1

Related Questions