# Matrix algebra
# rnorm function generates random variables from a normal distribution
rnorm(1)
# generate 10000 numbers and plot histogram
hist(rnorm(10000))

# same but for uniform distribution
hist(runif(10000))

# poisson distribution mean of 1
hist(rpois(10000, 1))

# make a vectors of random numbers
v1<-rnorm(10,mean = 3, sd = 1)
v2<-rnorm(10, mean = 20, sd= 5)
v3<-rnorm(10, mean = 1, sd = 0.5)
v4<-rnorm(10, mean = 30, sd= 1)
v5<-rnorm(10, mean = 2, sd= 5)

# combine rows into a matrix
rbind(v1,v2,v3,v4,v5)

# combine columns into a matrix
sample_matrix<-cbind(v1,v2,v3,v4,v5)

# produce a sequence of numbers
i<-seq(from=1,to=10)

# randomize the order of i
sample(i,10,replace=F)

# use that to randomize the order of the matrix
sample_matrix[sample(i,10,replace=F),]

# matrix functions
a<-matrix(c(1,4,-3,2,0,-1),ncol=3)
b<-matrix(c(6,7,3,-2,-3,5),ncol=3)
a+b
a*b

c<-matrix(c(1,2,3),ncol=1)
a %*% c

V1<-seq(1,6)
V2<-seq(1,10)
V1
V2
# none of these work...
V1+V2
V1*V2
V1 %*% V2

# matrix multiplication works and produces an outer product
V1 %*% t(V2)

# do not conform
t(V1) %*% V2

# inner product
t(V1) %*% V1

# recall our simple matrix from above
sample_matrix

# make a one-dimensional matrix of 1's as long as a column
one_vec<-rep(1,10)
one_vec

# column totals calculated through matrix multiplication 
t(one_vec) %*% sample_matrix

# divide by the number of samples to get column means
col_means<-t(one_vec) %*% sample_matrix/10

# should be some numbers as in colMeans function
colMeans(sample_matrix)

# center the matrix by column means
mean_centered<-sample_matrix - one_vec %*% col_means
mean_centered

# multiply the mean centered matrix by itself/6-1
varcovar<-t(mean_centered) %*% mean_centered/(10-1)
varcovar

# understand where these numbers come from...
# multiply row 1 in the first matrix by column one in the second matrix
t(mean_centered)[1,] * (mean_centered/(10-1))[,1]

# sum those numbers to get value [1,1] in the variance covariance matrix
sum(t(mean_centered)[1,] * (mean_centered/(10-1))[,1])

# should produce the same output as the cov function
cov(sample_matrix)

# calculating correlations - covariance/ square root of product variance
varcovar[1,2]/sqrt(varcovar[1,1]*varcovar[2,2])
# should be the same as this
cor(sample_matrix[,1],sample_matrix[,2])
# same as value [1,2] in a correlation matrix
cor(sample_matrix)