The code II

# create least squares system:
A = transpose(array([x, zeros(n, Float)+1]))
B = y
from LinearAlgebra import linear_least_squares
sol = linear_least_squares(A, B)
# sol is a 4-tuple, the solution (a,b) is the 1st entry:
a, b = sol[0]

import Gnuplot
g = Gnuplot.Gnuplot(persist=1)
# persist=1: let plot remain on the screen
g('set pointsize 2')
d1 = Gnuplot.Data(x, y,
     with='points', title='data')
d2 = Gnuplot.Func("%(a_exact)g*x + %(b_exact)g" % \
     vars(), with='lines', title='exact')
d3 = Gnuplot.Func("%(a)g*x + %(b)g" % vars(),
     with='lines', title='least-squares fit')
g.plot(d1, d2, d3)
g.hardcopy(filename='gp_tmp.ps', enhanced=1, 
  color=0, fontname='Times-Roman', fontsize=28)