using PyPlot

x = linspace(0, 3, 200)
y = sqrt(x.^2 + 1)                # hyperbola
yn = y + randn(length(y))*0.1     #  + noise
plot(x, y, "k--")
plot(x, yn, "r.")
legend(["hyperbola", "hyperbola + noise"], loc="upper left")

A = [x.^0 x.^1 x.^2] # .^n is element-wise exponentiation and [a b c] concatenates columns into a matrix

# More generally, .^ is a "broadcasting" operation that can combine a column vector x and a row vector [0 1 2]
# into a matrix:
A = x .^ [0 1 2]

using Interact

f = figure()
@manipulate for n = 1:40
    withfig(f) do
        A = x .^ [0:n]'
        c = A \ yn
        plot(x, yn, "r.")
        plot(x, A * c, "b-")
        legend(["data", "fit to degree $n"], loc="upper left")
    end
end

N = 1:20
semilogy(N, [cond(x .^ [0:n]') for n in N], "bo-")
xlabel(L"polynomial degree $n$")
ylabel(L"condition number of $A$")
title("condition number of least-square polynomial fitting")