topic Re: Help to solve equation in SAS Procedures

Help to solve equation

bollibompa — Thu, 25 Jun 2015 12:59:17 GMT

Hi,

I have issues writning a formula to solve the following equation in SAS.

I have the following 10 observations (below).

Var1 is unknown ( Var1 have the same for all observations)

Var2 is known.

The value for Var3 for each observation is calculated using the following formula (EXP(Var1+Var2)/(1+EXP(Var1+Var2))

Known is also that the sum of Var3 for all ten observations is 5. This means that Var1 is about -0.56.

I now want SAS to automatic solve X for me

Thanks for help!

/Thomas

obs	var1	var2	var3
1	X	0.1	EXP(X+0.1)/(1+EXP(X+0.1))
2	X	0.2	EXP(X+0.2)/(1+EXP(X+0.2))
3	X	0.3	EXP(X+0.3)/(1+EXP(X+0.3))
4	X	0.4	EXP(X+0.4)/(1+EXP(X+0.4))
5	X	0.5	EXP(X+0.5)/(1+EXP(X+0.5))
6	X	0.6	EXP(X+0.6)/(1+EXP(X+0.6))
7	X	0.7	EXP(X+0.7)/(1+EXP(X+0.7))
8	X	0.8	EXP(X+08)/(1+EXP(X+0.8))
9	X	0.9	EXP(X+0.9)/(1+EXP(X+0.9))
10	X	1.0	EXP(X+1.0)/(1+EXP(X+1.0))

Re: Help to solve equation

LinusH — Thu, 25 Jun 2015 13:02:30 GMT

You need to put you equation in your code, not as a variable. Should not be any problem since the equation logic is the same for all rows.

Re: Help to solve equation

Ksharp — Thu, 25 Jun 2015 14:03:28 GMT

It looks like that you are trying to look for a root .You might want post it at IML forum . Rick might have a good idea. Here is my GA code.

But My root X=-6.392141 not -0.56 .

Code: Program

data have;
input obs var1 $ var2 var3 : $40.;
cards;
1 X 0.1 EXP(X+0.1)/(1+EXP(X+0.1))
2 X 0.2 EXP(X+0.2)/(1+EXP(X+0.2))
3 X 0.3 EXP(X+0.3)/(1+EXP(X+0.3))
4 X 0.4 EXP(X+0.4)/(1+EXP(X+0.4))
5 X 0.5 EXP(X+0.5)/(1+EXP(X+0.5))
6 X 0.6 EXP(X+0.6)/(1+EXP(X+0.6))
7 X 0.7 EXP(X+0.7)/(1+EXP(X+0.7))
8 X 0.8 EXP(X+08)/(1+EXP(X+0.8))
9 X 0.9 EXP(X+0.9)/(1+EXP(X+0.9))
10 X 1.0 EXP(X+1.0)/(1+EXP(X+1.0))
;
run;

proc sql;
select var3 into : list separated by '+' from have;
quit;
%put &list ;


proc iml;
start func(x);
 v=abs(&list-5);
 return (v);
finish;

id=gasetup(1,1,1234);
call gasetobj(id,0,"func");
call gasetsel(id,100,1,0.95);
call gasetcro(id,0.05,4);
call gasetmut(id,0.05);

call gainit(id,10000,{-10,10});
niter = 100;
summary = j(niter,2);
mattrib summary [c = {"Min Value", "Avg Value"} l=""];
do i = 1 to niter;
 call garegen(id);
 call gagetval(value, id);
 summary[i,1] = value[1];
 summary[i,2] = value[:];
end;
call gagetmem(mem, value, id, 1);
print mem[ l="best member:"],
   "Min Value: " value[l = ""] ;
iteration = t(1:niter);
print iteration summary;
call gaend(id);
quit;

Xia Keshan

Re: Help to solve equation

bollibompa — Thu, 25 Jun 2015 14:22:08 GMT

Thanks for reply,

I have var3 in the code but I have problems to solve the equation due to the sum of all 10 observations:

Equation: sum of (EXP(Var1+Var2)/(1+EXP(Var1+Var2)) for all 10 observations =5

and I want SAS to solve Var1 for me (Var2 is known)

Thanks for further guidance if there are any procedures to solve equations loike this

/Thomas

Re: Help to solve equation

PGStats — Thu, 25 Jun 2015 19:50:33 GMT

It can be solved with proc fcmp :

proc fcmp outlib=sasuser.fcmp.test;

function myFunc(x, A

);

do i = 1 to dim(A);

sum + exp(x + A)/(1 + exp(x + A));

end;

return(sum);

endsub;

function solveForX(A

, sum, initX);

array solvopts[1] initial;

initial = initX;

return (solve("myFunc", solvopts, sum, ., A));

endsub;

run;

options cmplib=sasuser.fcmp;

data _null_;

array var2{10};

do i = 1 to 10;

var2{i} = i / 10;

output;

end;

expectedSum = 5;

initialGuess = -0.5;

x = solveForX(var2, expectedSum, initialGuess);

put x=;

run;

Re: Help to solve equation

Rick_SAS — Thu, 25 Jun 2015 20:08:10 GMT

This is a one-dimensional root-finding problem: Find the value of alpha such that

sum( logistic(alpha + var2) ) = 5

sum( logistic(alpha + var2) ) - 5 = 0.

Just do an internet search for "find root sas."

Re: Help to solve equation

bollibompa — Thu, 25 Jun 2015 22:30:50 GMT

Thanks thats brilliant!

A following question. In my example the values of var2 went from 0.1 to 1.

However, values for var2 might be different not following the a specific pattern (below). Is it possible to use the same proc fcmp-model and add the dataset I want to use to the array?

/Thomas

obs	var1	var2	var3
1	X	0.15	EXP(X+0.15)/(1+EXP(X+0.15))
2	X	0.03	EXP(X+0.03)/(1+EXP(X+0.03))
3	X	0.23	EXP(X+0.23)/(1+EXP(X+0.23))
4	X	0.44	EXP(X+0.44)/(1+EXP(X+0.44))
5	X	0.16	EXP(X+0.16)/(1+EXP(X+0.16))
6	X	0.7	EXP(X+0.7)/(1+EXP(X+0.7))
7	X	0.12	EXP(X+0.12)/(1+EXP(X+0.12))
8	X	0.32	EXP(X+0.32)/(1+EXP(X+0.32))
9	X	0.22	EXP(X+0.22)/(1+EXP(X+0.22))
10	X	0.82	EXP(X+0.82)/(1+EXP(X+0.82))

Re: Help to solve equation

PGStats — Fri, 26 Jun 2015 01:26:21 GMT

Replace the _null_ step with

/* Prepare equation parameters in a dataset */

data var2;

input eqNo var2_1-var2_10 expectedSum initialGuess;

datalines;

1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 5 -0.5

2 0.15 0.03 0.23 0.44 0.16 0.7 0.12 0.32 0.22 0.82 5 -0.5

;

/* solve equations */

data want;

set var2;

array var2Array{*} var2_:;

x = solveForX(var2Array, expectedSum, initialGuess);

run;

proc print data=want noobs; run;

Re: Help to solve equation

Ksharp — Fri, 26 Jun 2015 12:31:08 GMT

You mean creating a VAR3 variable ,right ?

Surprise! This time I got root X=-0.55 . There must be some wrong in your original var3.

Code: Program

data have;
input obs var1 $ var2;
var3=cats('EXP(X+',var2,')/(1+EXP(X+',var2,'))');
cards;
1 X 0.1 
2 X 0.2 
3 X 0.3 
4 X 0.4 
5 X 0.5 
6 X 0.6 
7 X 0.7 
8 X 0.8 
9 X 0.9 
10 X 1.0 
;
run;

proc sql;
select var3 into : list separated by '+' from have;
quit;
%put &list ;


proc iml;
start func(x);
 v=abs(&list-5);
 return (v);
finish;

id=gasetup(1,1,1234);
call gasetobj(id,0,"func");
call gasetsel(id,100,1,0.95);
call gasetcro(id,0.05,4);
call gasetmut(id,0.05);

call gainit(id,10000,{-10,10});
niter = 100;
summary = j(niter,2);
mattrib summary [c = {"Min Value", "Avg Value"} l=""];
do i = 1 to niter;
 call garegen(id);
 call gagetval(value, id);
 summary[i,1] = value[1];
 summary[i,2] = value[:];
end;
call gagetmem(mem, value, id, 1);
print mem[ l="best member:"],
   "Min Value: " value[l = ""] ;
iteration = t(1:niter);
print iteration summary;
call gaend(id);
quit;

Xia Keshan

Re: Help to solve equation

Ksharp — Fri, 26 Jun 2015 12:45:20 GMT

Here is How to find a root by Rick . But I don't agree with him a little bit .

http://blogs.sas.com/content/iml/2015/06/22/root-guess.html

Re: Help to solve equation

billfish — Fri, 26 Jun 2015 22:51:41 GMT

A pedestrian approach.

let us generalize a bit the function for which one is looking for roots.

The following applies 3 parameters (p1,p2,p3) to a slightly generalized function given in the Question.

************************************************;
**** if you want, create sets of parameters ****;
************************************************;
data t_param(keep=zSteps i p1 p2 p3);
   do k= 1 to 5;
      zSteps=10*(k=1)+20*(k=2)+40*(k=3)+50*(k=4)+100*(k=5);
   do i=1 to zSteps;
      p1=10/zSteps;
      p2=(1/zSteps)*i;
      p3=(1/zSteps)*i;
      output;
   end;
   end;
run;

**************************************************;
**** macro funk() can be an arbitrary function ***;
**************************************************;
%macro funk(d1,d2,d3);
    %global fnc;
    %let fnc = &d1*exp(x+&d2)/(1+exp(x+&d3));
    &fnc.
%mend;

*****************************************;
**** table t_root collects the roots ****;
*****************************************;
data t_root(keep=zSteps ySteps yRoot);
   array pp1(1000);
   array pp2(1000);
   array pp3(1000);
   array fnk(10000);
   ySteps=400;
   LHS_eqn=5; **** this is equation the Left-Hand-Side value ****;

   nn=0;
   call missing(of pp1(*));
   call missing(of pp2(*));
   call missing(of pp3(*));
   call missing(of fnk(*));
   *********************************************;
   **** load the matrices pp1() pp2() pp3() ****;
   *********************************************;
   do until (last.zSteps);
      nn+1;
      set t_param;
      by zSteps;
      pp1(nn)=p1; pp2(nn)=p2; pp3(nn)=p3;
   end;
   ********************************;
   **** find f(x) from -1 to 1 ****;
   ********************************;
   do yIter=1 to ySteps;
      x=-1.0+2*(yIter/ySteps);
      FF=-1*LHS_eqn;
   do c1=1 to nn;
      f1=pp1(c1); f2=pp2(c1); f3=pp3(c1);
      FF = FF + %funk(f1,f2,f3);
   end;
      fnk(yIter)=FF;
   end;
   *********************************;
   **** find roots from -1 to 1 ****;
   *********************************;
   do yIter=1 to (ySteps-1);
      sgn1=sign(fnk(yIter+1));
      sgn0=sign(fnk(yIter+0));
      if (sgn1=0) then do; yRoot=-1.0+2*((yIter+1)/ySteps); output; end;
      else if not(sgn0=0) and not(sgn1=sgn0) then do;
         yRoot= -1.0+2*((yIter+1)/ySteps);
         abs1=abs(fnk(yIter+1));
         abs0=abs(fnk(yIter+0));
         yRoot=yRoot -(2/ySteps)*(abs1/(abs0+abs1));
         output;
      end;
   end;
run;

table t_root becomes:

ySteps zSteps   yRoot
==========================
400      10     -0.5500
400      20     -0.5250
400      40     -0.5125
400      50     -0.5100
400     100     -0.5050
==========================

Re: Help to solve equation

bollibompa — Sat, 27 Jun 2015 06:51:43 GMT

Thank you all for your relies!

I tried the solution by PGstats and it worked well for me

thanks

thomas