topic how to get the point of intersection on two plots? in SAS Procedures

how to get the point of intersection on two plots?

TomiKong — Sat, 03 Mar 2012 09:57:07 GMT

For example, I used proc gplot to do plots on y vs x and z vs x. How to get get the point of intersection on two plots? Thanks.

data mydata;

input x y z;

datalines;

1 550 100

2 500 200

3 450 300

4 400 400

5 350 500

6 300 600

7 250 700

8 200 800

;

run;

symbol1 v=none INTERPOL=join c=blue width=5;

symbol2 v=none INTERPOL=join c=red width=5;

proc gplot data=mydata;

plot (y z)*x/overlay;

run;

quit;

how to get the point of intersection on two plots?

data_null__ — Sat, 03 Mar 2012 17:14:53 GMT

Your question seems somewhat vague. Clearly the lines intersect at x=4 y,z=400. Are you asking where the lines intersect if there is no common point in the data?

You would need to calculate that point, perhaps someone can show us the SAS Procedure for this.

how to get the point of intersection on two plots?

TomiKong — Sat, 03 Mar 2012 17:33:21 GMT

Thanks a lot for your reply. It's easy to get the intersection point in this example. But for some complicated plots, I need to find one way to calculate the point where plots intersect.

how to get the point of intersection on two plots?

art297 — Sat, 03 Mar 2012 18:18:56 GMT

I think they are addressing a similar problem in the following thread: http://www.listserv.uga.edu/cgi-bin/wa?A2=ind1002c&L=sas-l&D=1&O=A&J=0&P=31586

Hopefully, you can extract what you are looking for from at least one of the methods suggested.

Re: how to get the point of intersection on two plots?

data_null__ — Sat, 03 Mar 2012 18:34:51 GMT

I hate to post this program because I "know" there has to be an easy way to do this using a SAS PROC. Some of you statisticians should be able to help.

This programs plots two lines with the space between filled. Finding the points where the two lines cross is needed to make the fill work properly. The program uses JOIN interpolation, more complex interpolation (I don't know). You will recognize the formula for the intersection of two lines.

I used an associative array so I could do "everything" in one data step.

The Y values are transformed to get one Y*X set for the filled part of the graph.

This does what you want (find the intersection) and then applies it.

goptions reset=all;
data test;
   input x0 yA yB; * had to change your names a bit;
   datalines;
1          550          100
2          500          575
3          450          300
4          400          400
5          350          500
6          300          600
7          850          700
8          200          800
;
run;
 
data _null_;
   if 0 then do;
      set test;
      by x0;
      end;
   length difx dif1 dif2 xL y1L y2L m1 m2 b1 b2 _m1 _m2 _b1 _b2 x y f1-f3 8;
   length _name_ $32;
   dcl hash d(ordered:'Y');
   dcl hiter di('d');
   do while(1);
      call vnext(_name_);
      put _name_=;
      if upcase(_name_) eq '_NAME_' then leave;
      if _name_ eq: 'FIRST.' then do;
         d.definekey(scan(_name_,-1,'.'));
         continue;
         end;
      else if _name_ eq: 'LAST.' then continue;
      d.definedata(_name_);
      end;
   d.definedone();
   * load data into hash, insures ordered by x0;
   do until(eof);
      set test end=eof;
      rc = d.add();
      end;
   * iterate over the hash, as if it were an array;
   do rc = di.first() by 0 while(rc eq 0);
      *put _all_;
      difx = dif(x0);
      dif1 = dif(yA);
      dif2 = dif(yB);
 
      xl   = lag(x0);
      y1l  = lag(yA);
      y2l  = lag(yB);
 
      * compute slope intercept for each pair of points;
      m1   = divide(dif1,difx);
      b1 = yA - (m1*x0);
 
      m2   = divide(dif2,difx);
      b2 = yB - (m2*x0);
      
      * find the intersection of each pair of lines;
      *m1(x) + b1 = m2(x) + b2;
      _m1=m1; _b1=b1;
      _m2=m2; _b2=b2;
      * factor m2;
      select(sign(m2));
         when(-1) do; _m1=_m1+abs(_m2); _m2=0; end;
         when(1)  do; _m1=_m1-abs(_m2); _m2=0; end;
         otherwise _m2=0;
         end;
      * factor b1;
      select(sign(b1));
         when(-1) do; _b2=_b2+abs(_b1); _b1=0; end;
         when(1)  do; _b2=_b2-abs(_b1); _b1=0; end;
         otherwise _b1=0;
         end;
      * Intersection(x,y) of lines;
      x = divide(_b2,_m1);
      y = m1*x + b1;
 
      * do they cross within the bounds of the points.; 
      f1 = min(yA,y1L) lt y lt max(yA,y1L);
      f2 = min(yB,y2L) lt y lt max(yB,y2L);
      f3 = min(x0,xL)  lt x lt max(x0,xL); 
      f = d.replace();
      rc = di.next();
      end;
   * Output all these calculations as a check;
   d.output(dataset:'intersection');
   
   length type order 8;
   dcl hash area(ordered:'Y');
   area.definekey('type','order');
   area.definedata('type','order','x0','value','yA','yB');
   area.definedone();
 
   do type = 1;  
      do rc = di.first() by 0 while(rc eq 0);
         order = x0;
         value = max(yA,yB);
         rcA = area.add();
         if sum(f1,f2,f3) gt 0 then do;
            if round(x,1e-6) ne round(xL,1e-6) and round(x,1e-6) ne round(x0,1e-6)then do;
               x0   = x;
               order = x;
               value = y;
               call missing(yA,yB);
               rcA = area.add();
               end;
            end;
         rc = di.next();
         end;
      end;      
   do type = 2;  
      do rc = di.last() by 0 while(rc eq 0);
         order = -x0;
         value = min(yA,yB);
         call missing(yA,yB);
         rcA = area.add();
         if sum(f1,f2,f3) gt 0 then do;
            if round(x,1e-6) ne round(xL,1e-6) and round(x,1e-6) ne round(x0,1e-6)then do;
               order = -x;
               value = y;
               x0   = x;
               rcA = area.add();
               end;
            end;
         rc = di.prev();
         end;
      end;      
   rcA = area.output(dataset:'AREA');
   stop;
   run;
proc print data=intersection;
   run;
proc contents data=area;
proc print data=area;
   run;
symbol interpol=msolid cv=lightgreen co=darkgreen;
symbol2 v=dot  color=black i=join w=2;
symbol3 v=star color=red i=join w=2;
proc gplot data=area;
   plot value*x0=1 yA*x0=2 yB*x0=3 / overlay;
   run;
   quit;

Message was edited by: data _null_

Re: how to get the point of intersection on two plots?

PGStats — Sat, 03 Mar 2012 23:12:22 GMT

The interpolation specification I=JOIN generates line segments, so there could be many intersection between the two curves. You can find all the intersention points with the program (I added an x value to cover a special case) :

data mydata;
input x y z;
datalines;
0          600          600
1          550          100
2          500          200
3          450          300
4          400          400
5          350          500
6          300          600
7          250          700
8          200          800
;
run;

proc sort data=myData; by x; run;

data meet(keep=x1 x2 xi yzi);
set myData;
x1=lag(x); x2=x;
y1=lag(y); y2=y;
z1=lag(z); z2=z;
if _n_ > 1 then do;
alpha = (z2-y2) / (y1-y2-z1+z2);
if alpha ge 0 and alpha lt 1 then do;
    xi = alpha*x1 + (1-alpha)*x2;
    yzi = alpha*y1 + (1-alpha)*y2;
    output;
end;
end;
else if y=z then do;
xi = x;
yzi = y;
output;
end;
run;

/* at this point, intersection points coordinates are in dataset meet */

proc sql;
create table myGraphData as
select A.*, B.xi, B.yzi
from myData as A left join meet as B on A.x = B.x2;

proc sgplot data=myGraphData;
series x=x y=y;
series x=x y=z;
scatter x=xi y=yzi / markerattrs=(symbol=CircleFilled size=8);
run;

Extended by PG :

of if you want to use proc gplot :

proc sql;

create table myGData as

select * from

(select x, y, z, . as yzi from myData)

UNION

(select xi, ., ., yzi from meet)

order by yzi, x;

symbol1 v=none INTERPOL=join c=blue width=5;

symbol2 v=none INTERPOL=join c=red width=5;

symbol3 v=dot INTERPOL=none c=green width=5;

proc gplot data=myGData;

plot (y z yzi)*x/overlay;

run;

quit;

Re: how to get the point of intersection on two plots?

Ksharp — Mon, 05 Mar 2012 09:59:12 GMT

Just like PGStats . It is only suited for line.

data mydata;
input x y z;
datalines;
1          550          100
2          500          200
3          450          300
4          400          400
5          350          500
6          300          600
7          250          700
8          200          800
;
run;

data mean(keep=x mean_y);
 set mydata(obs=2) end=last;
 x=(z-y)*dif(x)/(dif(y)-dif(z))  + x ;
 mean_y=(z-y)*dif(y)/(dif(y)-dif(z))  + y ;
 if last then output;
run;
data want;
 set mydata mean;run;
symbol1 v=none INTERPOL=join c=blue width=5;
symbol2 v=none INTERPOL=join c=red width=5;
symbol3 v=dot INTERPOL=none c=green width=5;
proc gplot data=want;
plot (y z mean_y)*x/overlay;
run;
quit;

Ksharp