Contributor
Posts: 73

# Help to solve equation

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))
Super User
Posts: 5,876

## Re: Help to solve equation

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.

Data never sleeps
Contributor
Posts: 73

## Re: Help to solve equation

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

Super User
Posts: 10,766

## Re: Help to solve equation

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

Posts: 5,521

## Re: Help to solve equation

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;

PG

PG
Contributor
Posts: 73

## Re: Help to solve equation

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))
Posts: 5,521

## Re: Help to solve equation

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;

PG

PG
Contributor
Posts: 73

## Re: Help to solve equation

Thank you all for your relies!

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

thanks

thomas

Super User
Posts: 10,766

## Re: Help to solve equation

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

Super User
Posts: 10,766

## Re: Help to solve equation

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

SAS Super FREQ
Posts: 4,239

## Re: Help to solve equation

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

sum( logistic(alpha + var2) ) = 5

or

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

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

Contributor
Posts: 52

## Re: Help to solve equation

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
==========================

Discussion stats
• 11 replies
• 558 views
• 1 like
• 6 in conversation