with(stats): x := [1,2,3,4,5,6,7]: describe[mean](x): AVE := x -> stats[describe,mean](x): AVE(x); SD := x -> stats[describe,standarddeviation](x): SD(x); y := [6,7,5,4,3,1,2]: AVE(y); SD(y); r := (x,y) -> stats[describe,linearcorrelation](x,y): with(statplots): scatter := (x,y) -> stats[statplots,scatter2d](x,y): scatter(x,y); sdline := (x,d1,d2) -> AVE(d2)+sign(r(d1,d2))*(SD(d2)/SD(d1))*(x-AVE(d1)): with(plots): l1 := plot(sdline(u,x,y),u=0..10): scat := scatter(x,y): display({l1,scat}); yb := [2,1,4,3,7,5,6]: avex := AVE(x): sdx := SD(x): avey := AVE(yb): sdy := SD(yb): rb := r(x,yb): avex, sdx, avey, sdy, rb; evalf(rb); l2 := plot(sdline(u,x,yb),u=0..10): scat2 := scatter(x,yb): display({l2,scat2});