topic Re: find roots of a quadratic equation for each obs in SAS Programming

find roots of a quadratic equation for each obs

PamG — Sun, 05 Jan 2025 15:40:43 GMT

I have a data set with y and x. Is there a way to find the 2 roots of a quadratic equation for each obs?

The equation is y - 2*x*R + 3*x*R**2 =0.

Sample data is

data have;
input y x;
datalines;
1 2
2 4
5 6
8 9
;
run;

Re: find roots of a quadratic equation for each obs

Tom — Sun, 05 Jan 2025 15:53:14 GMT

The quadratic equation from Wikipedia:

So first convert your equation to the standard form.

Then just calculate the various parts. -b , 2a and sqrt(b**2-4ac).

Make sure to check if 4ac > b**2 because those cases have imaginary roots.

If it has real roots then build that by combining the parts twice, once using addition and once using subtraction.

Re: find roots of a quadratic equation for each obs

PamG — Sun, 05 Jan 2025 17:02:11 GMT

Thanks for your response. I was looking for an easier way to find the roots without having to write the equations in data step. For another project I will be having cubic polynomials. In that case how do I save the 3 roots in 3 variables? I saw some examples online who are using PROC MODEL or PROC IML. But they use some sort of starting numbers to find the roots.

Re: find roots of a quadratic equation for each obs

PaigeMiller — Sun, 05 Jan 2025 17:58:32 GMT

I was looking for an easier way to find the roots without having to write the equations in data step.

Hint: there are a gazillion different articles out on the internet regarding how to perform specific mathematical and statistical calculations, many of these articles written by @Rick_SAS . Searching for these are always a good place to start.

See Rick's 2011 article: https://blogs.sas.com/content/iml/2011/08/03/finding-the-root-of-a-univariate-function.html

Hint: in your second message, you said "I was looking for an easier way to find the roots without having to write the equations in data step. For another project I will be having cubic polynomials." Please tell us all of your needs in your first message from now on.

Re: find roots of a quadratic equation for each obs

sbxkoenk — Sun, 05 Jan 2025 21:10:26 GMT

Find the root of a function by using the SAS DATA step
By Rick Wicklin on The DO Loop August 19, 2020
https://blogs.sas.com/content/iml/2020/08/19/find-root-function-sas-data-step.html
A simple way to find the root of a function of one variable
By Rick Wicklin on The DO Loop February 5, 2014
https://blogs.sas.com/content/iml/2014/02/05/find-the-root-of-a-function.html
Finding the root of a univariate function
By Rick Wicklin on The DO Loop August 3, 2011
https://blogs.sas.com/content/iml/2011/08/03/finding-the-root-of-a-univariate-function.html
If you use SAS Viya and SAS Studio 5.2 or later AND your site licenses SAS/IML , then there's a code snippet called:
Find Roots of Nonlinear Equation : Enables you to find the roots of a function of one variable. Finding the root (or zero) of a function enables you to solve nonlinear equations.
Paper 3099-2019 Find a Function Root By Genetic Algorithm
Xia Ke Shan, iFRE Inc., Beijing, China
Kurt Bremser, Allianz Technology GmbH, Austria (Presenter)
https://support.sas.com/resources/papers/proceedings19/3099-2019.pdf
root finding, proc optmodel
https://communities.sas.com/t5/SAS-Procedures/root-finding-proc-optmodel/m-p/298259
Halley's method for finding roots
By Rick Wicklin on The DO Loop August 24, 2016
https://blogs.sas.com/content/iml/2016/08/24/halleys-method-finding-roots.html

Ciao, Koen

Re: find roots of a quadratic equation for each obs

PamG — Sun, 05 Jan 2025 22:45:45 GMT

Thanks for the link. I looked up and got some code. However I am getting only one root.

**Example 3: Solve for a: P=2*e1 + a + a*a;
data test;
input p e1;
datalines;
12 1
24 2
36 3
;
run;

/* FCMP routine to find function roots */
proc fcmp outlib=work.fcmp.test;
function fn(e1, a);
    return (2*e1 + a + a*a );
endsub;
function findZero(p, e1,init);
    array solvopts[1] initial (0);
    initial = init;
    return (solve("fn", solvopts, p, e1, .)); 
endsub;
run;

options cmplib=work.fcmp.test;
         
data roots;
set test;
format root 5.2;
root = findZero(p, e1, 0.1);
run;

Re: find roots of a quadratic equation for each obs

sbxkoenk — Sun, 05 Jan 2025 23:56:35 GMT

Sorry for still another link ... but I would do it this way :

Home > Analytics > SAS/IML> find roots of function

https://communities.sas.com/t5/SAS-IML-Software-and-Matrix/find-roots-of-function/td-p/794659

A quadratic function (equation) may have one, two, or zero roots.

If the discriminant is less than 0, then the quadratic has no real roots. (but the quadratic equation still has two complex-valued roots)
If the discriminant is equal to zero then the quadratic has equal roots. (two roots with the same value, called a double root)
If the discriminant is more than zero then it has 2 distinct roots.

Ciao,

Koen

Re: find roots of a quadratic equation for each obs

PaigeMiller — Mon, 06 Jan 2025 00:29:00 GMT

The link I gave does not have a solution that uses PROC FCMP. So, I don't know what code you are using or what link you got it from.

Re: find roots of a quadratic equation for each obs

Ksharp — Mon, 06 Jan 2025 03:01:46 GMT

Yes. Check @Rick_SAS blogs. Risk has written many blog about this topic.

For your this special topic, you could use PROC NLIN ,SAS/IML .SAS/OR to get root of function.

https://blogs.sas.com/content/iml/2022/05/04/bootstrap-nonlinear-regression.html

https://blogs.sas.com/content/iml/2015/06/08/fit-circle.html

https://blogs.sas.com/content/iml/2018/08/15/optimization-nonlinear-constraints.html

Here is some example:

First is from Rick

data have;
input y x;
datalines;
12 1
24 2
36 3
;
run;

 
proc nlin data=have plots(stats=none)=(fitplot); 
   parameters R=1;                
   model y = 2*x + R + R**2;
run;

Second is from Rick


proc iml;
start ObjF(R)  global(x, y);    /* x and y are fixed vectors of data */
   F = ssq(2*x + R + R**2 - y);            /* minimize F = SSQ differences */
   return(F);
finish;
use have;  read all var {x y};  close ;  /* read observed data */
R = {10};                           /* initial guess for (x0, y0) */
opt = {0 1};                           /* options=(min, print level) */
call nlptr(rc,result,"ObjF", r, opt);  /* find optimum (x0, y0) */
print result;
quit;

Third is from @RobPratt


proc optmodel;
   set OBS;
   num x {OBS};
   num y {OBS};
   read data have into OBS=[_N_] x y;
   var R >= 0;
   min F = sum {i in OBS} (2*x[i]+R+R^2  - y[i])^2;
   solve;
   print R;
quit;

Re: find roots of a quadratic equation for each obs

Ksharp — Mon, 06 Jan 2025 03:11:29 GMT

If R was a vector:


data have;
input y x;
datalines;
12 1
24 2
36 3
;
run;



proc optmodel;
   set OBS;
   num x {OBS};
   num y {OBS};
   read data have into OBS=[_N_] x y;
   var R{OBS} >= 0;
   min F = sum {i in OBS} (2*x[i]+R[i]+R[i]^2  - y[i])^2;
   solve ;

   print F R;
quit;

Re: find roots of a quadratic equation for each obs

Ksharp — Mon, 06 Jan 2025 06:56:16 GMT

Your object funtion
P=2*e1 + a + a*a;
is monotonic , so you can't expect to get two roots.

Re: find roots of a quadratic equation for each obs

Rick_SAS — Mon, 06 Jan 2025 11:52:17 GMT

Initially you said you wanted to solve

(3*x)*R**2 - (2*x)*R + y =0, which has polynomial coefficients {3*x, -2*x, y}

But later you use

1*R**2 + 1*R + 2*e1 = 0

Which is it?

In any case, to get the roots of a low-degree polynomial, use the POLYROOT function in SAS IML. Create a SAS data set that has the polynomial coefficients ordered from highest degree to lowest (the constant term). Read the coefficients into IML and loop over the coefficients.

A polynomial of degree d always has d complex roots. If you are interested only in the real roots, then find the roots for which the imaginary part is 0 (or very tiny).

The following code assumes that you want to solve the equation in your original post. It then creates a SAS data set that has ONLY the polynomial coefficients. The IML program extracts the REAL roots.

data have;
input y x;
datalines;
-1 2
-2 3
-5 6
-8 9
;

/* root of
   (3*x)*R**2 - 2*x*R  + y = 0 
   Let a = 3*x, b= -2*x, and c=y and use the quadratic equation
*/
data Coef;
keep c2 c1 c0;
set have;
c2 = 3*x;
c1 = -2*x;
c0 = y;
run;

proc iml;
delta = 1E-14;   /* small number to use as a custoff value */
use Coef;
read all var _num_ into coeff;  /* coefficients of R^2, R, constant */
close;

numPolys = nrow(coeff);
deg = ncol(coeff)-1;
roots = j(numPolys, deg, .);  /* poly of degree d has d complex roots */
do i = 1 to numPolys;
   r = polyroot( coeff[i,] );   /* returns all complex roots */
   print i r;
   /* store only the real roots for which r[,2]=0 */
   realIdx = loc( abs(r[,2]) < delta );
   numRealRts = ncol(realIdx);
   if numRealRts > 0 then 
      roots[i, 1:numRealRts] = rowvec( r[realIdx, 1] );  /* get the real parts */
end;
PolyNumber = T(1:numPolys);
RootNumber = 1:deg;
print PolyNumber roots[c=RootNumber];

You can use the same IML program to find roots of higher-degree polynomials. For example, if you want to find cubic roots, use the following DATA step to create the coefficients and use the same IML code:

data Coef;
input c3 c2 c1 c0;
datalines;
1 -1 2 1
1 -1 -4 4
1 1 -3 -3
;

Re: find roots of a quadratic equation for each obs

PamG — Mon, 06 Jan 2025 17:46:52 GMT

Thank you so much Rick!!! This works perfectly.