BookmarkSubscribeRSS Feed
TomiKong
Fluorite | Level 6

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;

6 REPLIES 6
data_null__
Jade | Level 19

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.  

TomiKong
Fluorite | Level 6

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.

art297
Opal | Level 21

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.

data_null__
Jade | Level 19

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(1do; _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(1do; _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_

PGStats
Opal | Level 21

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;

PG

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;

PG
Ksharp
Super User

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

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
  • 6 replies
  • 3093 views
  • 0 likes
  • 5 in conversation