pearson := proc( x, y )
   local n, X, Y, Sx, Sy, r, k;

   n := LinearAlgebra[Dimension](x);

   X, Y := 0, 0;
   for k from 1 to n do
      X := evalf(X + x[k]);
      Y := evalf(Y + y[k]);
   end do;

   X, Y := X/n, Y/n;

   Sx, Sy := 0, 0;
   for k from 1 to n do
      Sx := evalf(Sx + (x[k]-X)^2);
      Sy := evalf(Sy + (y[k]-Y)^2);
   end do;
   Sx, Sy := sqrt(Sx/n), sqrt(Sy/n);

   r := 0;
   for k from 1 to n do
      r := evalf( r + (x[k]-X)*(y[k]-Y) );
   end do;

   return r/((n-1)*Sx*Sy);

end proc;