Help using Base SAS procedures

Help to solve equation

Reply
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,256

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

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

Super User
Posts: 9,681

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;

x.png

Xia Keshan

Respected Advisor
Posts: 4,646

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

    Respected Advisor
    Posts: 4,646

    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: 9,681

    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;

    x.png

    Xia Keshan

    Super User
    Posts: 9,681

    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: 3,477

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

    Ask a Question
    Discussion stats
    • 11 replies
    • 505 views
    • 1 like
    • 6 in conversation