//　2種のデータから回帰直線を求め，交差する点を求める 2006/10/27 T.Ogawa　作成
//  回帰直線 y1(i) = a1 x(i) + b1
//  回帰直線 y1(i) = a2 x(i) + b2
//  の係数a1,a2,b1,b2を求める．

//         n
//  EE = Sum(y - y(i))^2
//       i = 1
//
//         n
//     = Sum(a*x(i)+b - y(i))^2
//       i = 1

//　よって EEの最小値を求めるために

//   d(EE)/da = 0
//   d(EE)/db = 0

//  を求めると, aa,bbは
//
//  [ aa     [ sum(x^2) sum(x)  ^(-1)  [ sum(xy)
//    bb ] =   sum(x)   n     ]      *   sum(y) ]
//
// として求められる．
//




//　入力データ
inputData = [ 1 3  5
              3 9  7
              5 25 9
              7 52 11 ];

x = inputData(:,1);
y1 = inputData(:,2);
y2 = inputData(:,3);

//係数aa,bbを求める．
 n = size(x,1);
 
 sum_x1y1 = 0;
 sum_y1 = 0;
 sum_x1y2 = 0;
 sum_y2 = 0;
 sum_x2 = 0;
 sum_x1 = 0;
 
 xone = ones(x);
 
  sum_x1 = xone'*x;
  sum_y1 = xone'*y1;
  sum_y2 = xone'*y2;
  sum_x2 = x'*x;
  sum_x1y1 = x'*y1;
  sum_x1y2 = x'*y2;
    
  Coef_Matrix1 = inv([sum_x2 sum_x1;sum_x1 n])*[sum_x1y1;sum_y1];
  Coef_Matrix2 = inv([sum_x2 sum_x1;sum_x1 n])*[sum_x1y2;sum_y2];
 
  aa1 = Coef_Matrix1(1);
  bb1 = Coef_Matrix1(2);
  aa2 = Coef_Matrix2(1);
  bb2 = Coef_Matrix2(2);
  
// グラフ表示
plot(x,y1,'O',x,y2,'x',x,aa1*x+bb1,x,aa2*x+bb2)

f=get("current_figure");
a=f.children; // the handle on the Axes child
a.data_bounds=[min(x)-abs(min(x))*0.1 max(x)+abs(max(x))*0.1 min([y1;y2])-abs(min([y1;y2]))*0.1 max([y1;y2])+abs(max([y1;y2]))*0.1]; // xaxis_left xaxis_right yaxis_down yaxis_up 
a.x_label.text = "x";
a.y_label.text = "y";
xgrid(2)

if abs(aa1-aa2)>0,
 x_cross = (bb2-bb1)/(aa1-aa2)
 y_cross = aa1*x_cross+bb1
end