# Time-stamp: <2012-11-28 15:16 Petri Koistinen>

# Tämä on R-koodia

# Simuloidaan vektoreita 2-ulotteisesta normaalijakaumasta

# Odotusarvovektori
mu <- matrix(c(1, 3), 2, 1)

# Kovarianssimatriisin määräävät seuraavat parametrit
rho <- 0.3 # luku väliltä (0, 1)
sx <- 1.2 # sx > 0
sy <- 3.4 # sy > 0

# Muodostetaan kovarianssimatriisi
S <- matrix(c(sx^2, rho*sx*sy, rho*sx*sy, sy^2), 2, 2)

# Anetaan R:n laskea Choleskyn hajotelma, vaikka tunnemmekin
# analyyttiset kaavat

A <- t(chol(S))

# ... tarvittiin transponointi. Nyt S = A A^T.

print(S)
print(A %*% t(A))

# Simuloidaan yksi arvo

u <- matrix(rnorm(2), 2, 1)

print(x <- A %*% u + mu)


# Yhdestä arvosta ei ole hyötyä, jos jakaumaa pitäisi visualisoida.
# Simuloidaan N arvoa:

N <- 1000000

umat <- matrix(rnorm(2 * N), 2, N)
xmat <- A %*% umat + matrix(rep(mu, N), 2, N)

# Estimoidaan jakauman tunnusluvut

rowMeans(xmat)
cov(t(xmat))

# Visualisoidaan n ensimmäistä

n <- 2000

plot(xmat[1, 1:n], xmat[2, 1:n], xlab = 'x', ylab = 'y', pch = '.')