using PyPlot

x = randn(100) # 100 gaussian random numbers:
plot(x, zeros(x), "r.")
grid()

x = randn(10000)
plt[:hist](x, bins=100);

mean(x)

sum((x.-mean(x)).^2)/(length(x)-1)

var(x)

xbig = randn(10^7)
mean(xbig), var(xbig)

x = sin.(linspace(0,4π,100) .+ randn(100)*0.1) .* (1 .+ 0.3*randn(100))
y = sin.(linspace(0,4π,100) .+ randn(100)*0.1) .* (1 .+ 0.3*randn(100))
z = randn(100)
plot(x, "b.-")
plot(y, "r.-")
plot(z, "g.-")
legend(["x","y","z"])

mean(x),mean(y),mean(z)

# A simple covariance function.  See https://github.com/JuliaStats/StatsBase.jl for
# better statistical functions in Julia.
covar(x,y) = dot(x .- mean(x), y .- mean(y)) / (length(x) - 1)

covar(x,y), covar(x,z)

covar(x,y) / sqrt(covar(x,x) * covar(y,y)) # correlation, manually computed

cor(x,y)

cor(x,z)

abs(cor(x,y)/cor(x,z))

plot(x,y, ".")
xlabel("x")
ylabel("y")
axhline(0, linestyle="--", color="k", alpha=0.2)
axvline(0, linestyle="--", color="k", alpha=0.2)

A = [x' .- mean(x); y' .- mean(y)] # rows are x and y with means subtracted

S = A * A' / (length(x)-1)

σ², Q = eig(S)

σ²

plot(x.-mean(x),y.-mean(y), ".")
xlabel("x")
ylabel("y")
axhline(0, linestyle="--", color="k", alpha=0.2)
axvline(0, linestyle="--", color="k", alpha=0.2)
arrow(0,0, Q[:,1]..., head_width=0.1, color="r")
arrow(0,0, Q[:,2]..., head_width=0.1, color="r")
text(Q[1,1]+0.1,Q[2,1], "q₁", color="r")
text(Q[1,2]-0.2,Q[2,2], "q₂", color="r")
axis("equal")

U, σ, V = svd(A / sqrt(length(x)-1))

σ.^2

U

plot(x.-mean(x),y.-mean(y), ".")
xlabel("x")
ylabel("y")
axhline(0, linestyle="--", color="k", alpha=0.2)
axvline(0, linestyle="--", color="k", alpha=0.2)
arrow(0,0, U[:,1]..., head_width=0.1, color="r")
arrow(0,0, U[:,2]..., head_width=0.1, color="r")
text(U[1,2]+0.1,U[2,1], "u₁", color="r")
text(U[1,2]+0.2,U[2,2], "u₂", color="r")
axis("equal")