Estimation of nonnegative parameters

[Marx]

Dear Xia Keshan,

If X1, X2, X2, X4 <1, then I would like to know how to add these constraints in your coding. 

Thank you very much in advance.

Sincerely yours,

J1

[Ksharp]

data have;
infile cards truncover expandtabs;
input C	M	G	T;
cards;
2.36361198	1.69019608	2.549003262	1990
2.390935107	1.770852012	2.551449998	1991
2.431363764	1.84509804	2.551449998	1992
2.436162647	1.875061263	2.552668216	1993
2.530199698	2.064457989	2.555094449	1994
2.643452676	2.120573931	2.557507202	1995
2.722633923	2.1430148	2.559906625	1996
2.764176132	2.152288344	2.561101384	1997
2.794488047	2.146128036	2.564666064	1998
2.827369273	2.220108088	2.565847819	1999
2.870988814	2.352182518	2.568201724	2000
2.907948522	2.387389826	2.56937391	2001
2.938519725	2.469822016	2.571708832	2002
2.972665592	2.615950052	2.575187845	2003
3.024485668	2.748962861	2.57634135	2004
3.083860801	2.820857989	2.579783597	2005
3.151982395	2.900367129	2.582063363	2006
3.240049772	2.983626287	2.584331224	2007
3.343999069	3.058426024	2.586587305	2008
3.394101302	3.001733713	2.587710965	2009
3.457427693	3.145817714	2.591064607	2010
;
run;

proc iml;
use have;
read all var {C M G T};
close; 

start func(x) global(C,M,G,T);
 obj=sum(  (1/x[1])# ( (C#(M##x[2])# ((380/G)##x[3]) )##x[1] )  #exp(-x[4]#T)  );
 return (obj);
finish;
start C_UC2D(x);
 c = j(6,1,0); /*<-- 6 constraint*/
 c[1]=1-x[1];
 c[2]=1-x[2];
 c[3]=1-x[3];
 c[4]=1-x[4];         /*<-- Change it for new constraint*/
 c[5]=1-X[1]*(X[2]+X[3]);
 c[6]=1-X[1]*(1+X[2]+X[3]);
return(c);
finish;

call randseed(123456);
x = j(10,4,.); /*produce 10 solutions*/
temp=j(1,4,.); /*hold a temporary solution*/

opt = j(1, 11, .);
opt[1]=1;opt[2] = 0; opt[10] = 6;  /*<-- 6 constraint*/



do i=1 to nrow(x);
 call randgen(temp,'uniform');
 CALL NLPNMS(rc,xres,"func",temp,opt, , , , , "C_UC2D");
 x[i,]=xres; /*Save a solution*/
end;

create want from x[c=('x1':'x4')]; /*Save the solution into a table*/
append from x;
close;

quit;

proc print noobs;run;

[Marx]

Dear Xia Keshan,

Thank you very much for your response. 

To get standard errors, I added these code in your coding. 

proc means data=want mean stderr nway;
class parameter ;
var estimate;
run;

But results are 

SAS Output

Analysis Variable : Estimate Parameter N Obs Mean Std Error X1 14 X2 14 X3 14 X4 14

0.4280902	0.1635312
3.1404260	0.7945652
3.5133741	0.9782815
4.5563504	0.8566756

And I would like to get X2, X3, X4 are too big. I would like to know how to get standard errors in your coding. 

Thank you very much in advance. 

Sincerely your,

J1

[Ksharp]

Table WANT only contain 10 solutions. Try this code to get them std err:


proc means data=want mean stderr ;
run;

[Marx]

Dear Xia Keshan,

Thank you very much for your coding. 

I have one more question about your coding. 

If I want to 100 solutions, then I would like to know how to do it. 

I tried to change 10 to 100 below. But the result shows just empty box.

call randseed(123456);
x = j(100,4,.); /*produce 10 solutions*/
temp=j(1,4,.); /*hold a temporary solution*/

opt = j(1, 11, .);
opt[1]=1;opt[2] = 0; opt[100] = 6;  /*<-- 6 constraint*/

Sincerely yours,

J1

[Marx]

Dear Xia Keshan,

I ran this coding, and I found that the sixth constraint does not work (c[6]=1-X[1]*(1+X[2]+X[3])). i.e. X[1]*(1+X[2]+X[3]) = 3.72767...

"X[1]*(1+X[2]+X[3])" must be less than 1.

I would like to know hwo to fix this problem.  

Also, I would like to know if I add "bounds" condition in SAS/IML such as NLIN or not. 

Thank you very much in advance. 

Sincerley yours,

J1

data have;
infile cards truncover expandtabs;
input C M G T;
cards;
231 49 1.073446328 1990
246 59 1.06741573 1991
270 70 1.06741573 1992
273 75 1.06442577 1993
339 116 1.058495822 1994
440 132 1.052631579 1995
528 139 1.046831956 1996
581 142 1.043956044 1997
623 140 1.035422343 1998
672 166 1.032608696 1999
743 225 1.027027027 2000
809 244 1.02425876 2001
868 295 1.018766756 2002
939 413 1.010638298 2003
1,058 561 1.00795756 2004
1,213 662 1 2005
1,419 795 0.994764398 2006
1,738 963 0.989583333 2007
2,208 1,144 0.984455959 2008
2,478 1,004 0.981912145 2009
2,867 1,399 0.974358974 2010
;
run;

proc iml;
use have;
read all var {C M G T};
close;

start func(x) global(C,M,G,T);
obj=sum( (1/x[1])# ( (C#(M##x[2])# ((380/G)##x[3]) )##x[1] ) #exp(-x[4]#T) );
return (obj);
finish;
start C_UC2D(x);
c = j(6,1,0); /*<-- 6 constraint*/
c[1]=1-x[1];
c[2]=1-x[2];
c[3]=1-x[3];
c[4]=1-x[4]; /*<-- Change it for new constraint*/
c[5]=1-X[1]*(X[2]+X[3]);
c[6]=1-X[1]*(1+X[2]+X[3]);
return(c);
finish;

call randseed(123456);
x = j(10,4,.); /*produce 10 solutions*/
temp=j(1,4,.); /*hold a temporary solution*/

opt = j(1, 11, .);
opt[1]=1;opt[2] = 0; opt[10] = 6; /*<-- 6 constraint*/

do i=1 to nrow(x);
call randgen(temp,'uniform');
CALL NLPNMS(rc,xres,"func",temp,opt, , , , , "C_UC2D");
x[i,]=xres; /*Save a solution*/
end;

create want from x[c=('x1':'x4')]; /*Save the solution into a table*/
append from x;
close;

quit;

proc print noobs;run;

proc means data=want mean stderr ;
run;

[Ksharp]

Due to your ill Object function or Constrain condition, When I use 
x = j(100,4,.); /*produce 100 solutions*/


I will get these error:

 WARNING: The point x is feasible only at the LCEPSILON= 1E-7 range.
 NOTE: ABSXCONV convergence criterion satisfied.
 NOTE: ABSXCONV convergence criterion satisfied.
 NOTE: ABSCONV convergence criterion satisfied.
 NOTE: ABSXCONV convergence criterion satisfied.
 NOTE: ABSXCONV convergence criterion satisfied.
 ERROR: Overflow error in NLPNMS.
 

So there are must something wrong out there.
If you want check constrain condition. Just get a solution one time.
x = j(1,4,.); /*produce 100 solutions*/
I get these:
x1	x2	x3	x4
0.59408	0.51547	0.16519	0.72280

In order to check the sixth constraint condition,I wrote this:
proc iml;
start C_UC2D(x);
 c=X[1]*(1+X[2]+X[3]);
return(c);
finish;

x={0.59408,0.51547,0.16519,0.72280};
obj=C_UC2D(x);
print obj;
quit;

And got this :
obj
0.9984465

So it should be all right .

[Ksharp]

If you want check another solution ,you can change SEED:
call randseed(123456);
-->
call randseed(1);


And still get one solution one time:
do i=1 to nrow(x);





CODE:



data have;
infile cards truncover expandtabs;
input C	M	G	T;
cards;
2.36361198	1.69019608	2.549003262	1990
2.390935107	1.770852012	2.551449998	1991
2.431363764	1.84509804	2.551449998	1992
2.436162647	1.875061263	2.552668216	1993
2.530199698	2.064457989	2.555094449	1994
2.643452676	2.120573931	2.557507202	1995
2.722633923	2.1430148	2.559906625	1996
2.764176132	2.152288344	2.561101384	1997
2.794488047	2.146128036	2.564666064	1998
2.827369273	2.220108088	2.565847819	1999
2.870988814	2.352182518	2.568201724	2000
2.907948522	2.387389826	2.56937391	2001
2.938519725	2.469822016	2.571708832	2002
2.972665592	2.615950052	2.575187845	2003
3.024485668	2.748962861	2.57634135	2004
3.083860801	2.820857989	2.579783597	2005
3.151982395	2.900367129	2.582063363	2006
3.240049772	2.983626287	2.584331224	2007
3.343999069	3.058426024	2.586587305	2008
3.394101302	3.001733713	2.587710965	2009
3.457427693	3.145817714	2.591064607	2010
;
run;

proc iml;
use have;
read all var {C M G T};
close; 

start func(x) global(C,M,G,T);
 obj=sum(  (1/x[1])# ( (C#(M##x[2])# ((380/G)##x[3]) )##x[1] )  #exp(-x[4]#T)  );
 return (obj);
finish;
start C_UC2D(x);
 c = j(6,1,0); /*<-- 6 constraint*/
 c[1]=1-x[1];
 c[2]=1-x[2];
 c[3]=1-x[3];
 c[4]=1-x[4];         /*<-- Change it for new constraint*/
 c[5]=1-X[1]*(X[2]+X[3]);
 c[6]=1-X[1]*(1+X[2]+X[3]);
return(c);
finish;

call randseed(1);
x = j(1,4,.); /*produce 1 solutions*/
temp=j(1,4,.); /*hold a temporary solution*/

opt = j(1, 11, .);
opt[1]=1;opt[2] = 0; opt[10] = 6;  /*<-- 6 constraint*/



do i=1 to nrow(x);
 call randgen(temp,'uniform');
 CALL NLPNMS(rc,xres,"func",temp,opt, , , , , "C_UC2D");
 x[i,]=xres; /*Save a solution*/
end;

create want from x[c=('x1':'x4')]; /*Save the solution into a table*/
append from x;
close;

quit;

proc print noobs;run;

[Marx]

Dear Xia Keshan,

Can you run with these data below?

I still found that the constrains do not work. 

231 49 1.073446328 1990
246 59 1.06741573 1991
270 70 1.06741573 1992
273 75 1.06442577 1993
339 116 1.058495822 1994
440 132 1.052631579 1995
528 139 1.046831956 1996
581 142 1.043956044 1997
623 140 1.035422343 1998
672 166 1.032608696 1999
743 225 1.027027027 2000
809 244 1.02425876 2001
868 295 1.018766756 2002
939 413 1.010638298 2003
1,058 561 1.00795756 2004
1,213 662 1 2005
1,419 795 0.994764398 2006
1,738 963 0.989583333 2007
2,208 1,144 0.984455959 2008
2,478 1,004 0.981912145 2009
2,867 1,399 0.974358974 2010

Sincerely yours,

J1

[Ksharp]

OK. In order to test the sixth constraint condition, I change some code:




data have;
infile cards truncover expandtabs;
input C	M	G	T;
cards;
231 49 1.073446328 1990
246 59 1.06741573 1991
270 70 1.06741573 1992
273 75 1.06442577 1993
339 116 1.058495822 1994
440 132 1.052631579 1995
528 139 1.046831956 1996
581 142 1.043956044 1997
623 140 1.035422343 1998
672 166 1.032608696 1999
743 225 1.027027027 2000
809 244 1.02425876 2001
868 295 1.018766756 2002
939 413 1.010638298 2003
1,058 561 1.00795756 2004
1,213 662 1 2005
1,419 795 0.994764398 2006
1,738 963 0.989583333 2007
2,208 1,144 0.984455959 2008
2,478 1,004 0.981912145 2009
2,867 1,399 0.974358974 2010
;
run;

proc iml;
use have;
read all var {C M G T};
close; 

start func(x) global(C,M,G,T);
 obj=sum(  (1/x[1])# ( (C#(M##x[2])# ((380/G)##x[3]) )##x[1] )  #exp(-x[4]#T)  );
 return (obj);
finish;
start C_UC2D(x);
 c = j(6,1,0); /*<-- 6 constraint*/
 c[1]=1-x[1];
 c[2]=1-x[2];
 c[3]=1-x[3];
 c[4]=1-x[4];         /*<-- Change it for new constraint*/
 c[5]=1-X[1]*(X[2]+X[3]);
 c[6]=1-X[1]*(1+X[2]+X[3]);
return(c);
finish;

call randseed(123456);
x = j(10,4,.); /*produce 10 solutions*/
temp=j(1,4,.); /*hold a temporary solution*/
return_code=1:9;

opt = j(1, 11, .);
opt[1]=1;opt[2] = 0; opt[10] = 6;  /*<-- 6 constraint*/



do i=1 to nrow(x);
 call randgen(temp,'uniform');
 CALL NLPNMS(rc,xres,"func",temp,opt, , , , , "C_UC2D");
 if any(return_code=rc) then x[i,]=xres; /*Save a solution*/
end;


/****Check the sixth constraint condition***/
start check_six(x);
 c=X[1]*(1+X[2]+X[3]);
return(c);
finish;
six=j(nrow(x),5,.);
do i=1 to nrow(x);
 temp=x[i,];
 obj=check_six(temp);
 six[i,]=temp||obj;
end;
print six[c=(('x1':'x4')||'obj')];
/**********/



create want from x[c=('x1':'x4')]; /*Save the solution into a table*/
append from x;
close;

quit;




Here is what I got:
six
x1	x2	x3	x4	obj
0.5940822	0.5154669	0.1651917	0.722802	0.9984493
0.5496765	0.7738273	0.0372558	0.8633171	0.9955098
0.4999731	0.3276918	0.1576245	-0.234566	0.7426183
0.4874287	0.2300376	0.4407849	0.8105903	0.8144068
0.5630604	0.092952	0.5302821	0.854167	0.9139788
0.2705592	1	0.9719575	-0.236711	0.8040905
0.4986595	0.8796288	0.1257641	0.4467502	1.0000081
0.367851	0.5049997	0.8690966	0.5191776	0.8733137
0.7423242	0.1527639	0.1943564	0.5509057	1
0.1890787	0.0202677	0.9160893	0.5739596	0.3661238


The last column OBJ is the value of the sixth constraint condition. 
No problem at all.
Post that problem solution, if you think which one of solution is not right .