Estimation of nonnegative parameters

[Marx]

Dear Xia Keshan,

Thank you very much for your help.

I ran two cases (US and China) below with your coding, but I got the same results with the same standard errors.

There are two different data (US and China). I would like to know if this is correct or not. 

US case;

data have;
infile cards truncover expandtabs;
input C	M	G	T;
cards;
4,774 630 1.073446328 1990
4,964 624 1.06741573 1991
5,265 668 1.06741573 1992
5,545 720 1.06442577 1993
5,851 813 1.058495822 1994
6,129 903 1.052631579 1995
6,445 964 1.046831956 1996
6,785 1,060 1.043956044 1997
7,175 1,120 1.035422343 1998
7,665 1,250 1.032608696 1999
8,237 1,480 1.027027027 2000
8,648 1,400 1.02425876 2001
9,035 1,430 1.018766756 2002
9,521 1,550 1.010638298 2003
10,129 1,800 1.00795756 2004
10,774 2,030 1 2005
11,394 2,240 0.994764398 2006
11,960 2,370 0.989583333 2007
12,382 2,560 0.984455959 2008
12,289 1,980 0.981912145 2009
12,724 2,360 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;

proc means data=want mean stderr ;
run;

China case;

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;

proc means data=want mean stderr ;
run;

[Ksharp]

You miss the code which creating a table. And I got different Std Err.
If you want get different solution, Plz change SEED:
call randseed(123456);
call randseed(56);
call randseed(6);








data have;
infile cards truncover expandtabs;
input C	M	G	T;
cards;
4,774 630 1.073446328 1990
4,964 624 1.06741573 1991
5,265 668 1.06741573 1992
5,545 720 1.06442577 1993
5,851 813 1.058495822 1994
6,129 903 1.052631579 1995
6,445 964 1.046831956 1996
6,785 1,060 1.043956044 1997
7,175 1,120 1.035422343 1998
7,665 1,250 1.032608696 1999
8,237 1,480 1.027027027 2000
8,648 1,400 1.02425876 2001
9,035 1,430 1.018766756 2002
9,521 1,550 1.010638298 2003
10,129 1,800 1.00795756 2004
10,774 2,030 1 2005
11,394 2,240 0.994764398 2006
11,960 2,370 0.989583333 2007
12,382 2,560 0.984455959 2008
12,289 1,980 0.981912145 2009
12,724 2,360 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;

create want from x[c=('x1':'x4')];
append from x;
close;
quit;

title 'US';
proc print data=want noobs;run;
proc means data=want mean stderr ;
run;






/*************************/
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;

create want from x[c=('x1':'x4')];
append from x;
close;
quit;

title 'China';
proc print data=want noobs;run;
proc means data=want mean stderr ;
run;





OUTPUT:

US
x1	x2	x3	x4
0.59408	0.51547	0.16519	0.72280
0.54968	0.77383	0.03726	0.86332
0.44058	0.32769	0.15762	0.19401
0.48743	0.23004	0.44078	0.81059
0.56306	0.09295	0.53028	0.85417
0.22858	0.76381	0.60680	0.02690
0.49866	0.87963	0.12576	0.44675
0.36785	0.50500	0.86910	0.51918
0.74232	0.15276	0.19436	0.55091
0.18908	0.02027	0.91609	0.57396
US
The MEANS Procedure
Variable	Mean	Std Error
x1
x2
x3
x4
0.4661315
0.4261442
0.4043245
0.5562577
0.0531292
0.0973758
0.1005858
0.0882303
China
x1	x2	x3	x4
0.59408	0.51547	0.16519	0.72280
0.54968	0.77383	0.03726	0.86332
0.49997	0.32769	0.15762	-0.23457
0.48743	0.23004	0.44078	0.81059
0.56306	0.09295	0.53028	0.85417
0.27056	1.00000	0.97196	-0.23671
0.49866	0.87963	0.12576	0.44675
0.36785	0.50500	0.86910	0.51918
0.74232	0.15276	0.19436	0.55091
0.18908	0.02027	0.91609	0.57396
China
The MEANS Procedure
Variable	Mean	Std Error
x1
x2
x3
x4
0.4762694
0.4497636
0.4408403
0.4870393
0.0511342
0.1086841
0.1144286
0.1289037

[Ksharp]

OK. You are right. After adjusting the data step to input data correctly . I got the same output.
What I can say is both solution should be right. 
If you want get different solution ,change SEED as I said before.





data have;
infile cards truncover expandtabs;
input (C 	M	G	T) (: comma32.);
cards;
4,774 630 1.073446328 1990
4,964 624 1.06741573 1991
5,265 668 1.06741573 1992
5,545 720 1.06442577 1993
5,851 813 1.058495822 1994
6,129 903 1.052631579 1995
6,445 964 1.046831956 1996
6,785 1,060 1.043956044 1997
7,175 1,120 1.035422343 1998
7,665 1,250 1.032608696 1999
8,237 1,480 1.027027027 2000
8,648 1,400 1.02425876 2001
9,035 1,430 1.018766756 2002
9,521 1,550 1.010638298 2003
10,129 1,800 1.00795756 2004
10,774 2,030 1 2005
11,394 2,240 0.994764398 2006
11,960 2,370 0.989583333 2007
12,382 2,560 0.984455959 2008
12,289 1,980 0.981912145 2009
12,724 2,360 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(14,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;

create want from x[c=('x1':'x4')];
append from x;
close;
quit;

title 'US';
proc print data=want noobs;run;
proc means data=want mean stderr ;
run;






/************************/
data have;
infile cards truncover expandtabs;
input (C	M	G	T) (: comma32.);
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);
 c[1]=1-x[1];
 c[2]=1-x[2];
 c[3]=1-x[3];
 c[4]=1-x[4];      
 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(14,4,.);
temp=j(1,4,.);
return_code=1:9;

opt = j(1, 11, .);
opt[1]=1;opt[2] = 0; opt[10] = 6;  



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; 
end;

create want from x[c=('x1':'x4')];
append from x;
close;
quit;

title 'China';
proc print data=want noobs;run;
proc means data=want mean stderr ;
run;

[Ksharp]

After changing SEED , I got different Std Err.
But it is not trust due to the fact there are every small (10) number of solutions.
You need to get hundreds and thousands of solution to estimate std err.
But your ill object function or constraint condition don't allow it happend.








data have;
infile cards truncover expandtabs;
input (C 	M	G	T) (: comma32.);
cards;
4,774 630 1.073446328 1990
4,964 624 1.06741573 1991
5,265 668 1.06741573 1992
5,545 720 1.06442577 1993
5,851 813 1.058495822 1994
6,129 903 1.052631579 1995
6,445 964 1.046831956 1996
6,785 1,060 1.043956044 1997
7,175 1,120 1.035422343 1998
7,665 1,250 1.032608696 1999
8,237 1,480 1.027027027 2000
8,648 1,400 1.02425876 2001
9,035 1,430 1.018766756 2002
9,521 1,550 1.010638298 2003
10,129 1,800 1.00795756 2004
10,774 2,030 1 2005
11,394 2,240 0.994764398 2006
11,960 2,370 0.989583333 2007
12,382 2,560 0.984455959 2008
12,289 1,980 0.981912145 2009
12,724 2,360 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(654321);
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;

create want from x[c=('x1':'x4')];
append from x;
close;
quit;

title 'US';
proc print data=want noobs;run;
proc means data=want mean stderr ;
run;






/************************/
data have;
infile cards truncover expandtabs;
input (C	M	G	T) (: comma32.);
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);
 c[1]=1-x[1];
 c[2]=1-x[2];
 c[3]=1-x[3];
 c[4]=1-x[4];      
 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,.);
temp=j(1,4,.);
return_code=1:9;

opt = j(1, 11, .);
opt[1]=1;opt[2] = 0; opt[10] = 6;  



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; 
end;

create want from x[c=('x1':'x4')];
append from x;
close;
quit;

title 'China';
proc print data=want noobs;run;
proc means data=want mean stderr ;
run;