BookmarkSubscribeRSS Feed
bollibompa
Quartz | Level 8

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

11 REPLIES 11
LinusH
Tourmaline | Level 20

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
bollibompa
Quartz | Level 8

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

Ksharp
Super User

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

PGStats
Opal | Level 21

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
    bollibompa
    Quartz | Level 8

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

    PGStats
    Opal | Level 21

    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
    bollibompa
    Quartz | Level 8

    Thank you all for your relies!

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

    thanks

    thomas

    Ksharp
    Super User

    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

    Ksharp
    Super User

    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

    Rick_SAS
    SAS Super FREQ

    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."

    billfish
    Quartz | Level 8

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

    sas-innovate-2024.png

    Join us for SAS Innovate April 16-19 at the Aria in Las Vegas. Bring the team and save big with our group pricing for a limited time only.

    Pre-conference courses and tutorials are filling up fast and are always a sellout. Register today to reserve your seat.

     

    Register now!

    What is Bayesian Analysis?

    Learn the difference between classical and Bayesian statistical approaches and see a few PROC examples to perform Bayesian analysis in this video.

    Find more tutorials on the SAS Users YouTube channel.

    Click image to register for webinarClick image to register for webinar

    Classroom Training Available!

    Select SAS Training centers are offering in-person courses. View upcoming courses for:

    View all other training opportunities.

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