R9. Correlation in R

x and y correspond to 50 thousand samples of 2 throws of dice.

The correlation is near 0, as the two samples are independent, and will go to 0 as the number of samples is increased.

The calculated correlation and the result from cor(x,y) are compared.

# ex9.R
num_samples <- 50000
x <- sample(6, num_samples, replace = T)
y <- sample(6, num_samples, replace = T)
df <- data.frame(x,y)
df["a"] <- df["x"]-mean(df[["x"]])
df["b"] <- df["y"]-mean(df[["y"]])
df["ab"] <- df["a"]*df["b"]
df["sqa"] <- df["a"]*df["a"]
df["sqb"] <- df["b"]*df["b"]
den <- sqrt(sum(df[["sqa"]])*sum(df[["sqb"]]))
correlation <- sum(df[["ab"]])/den
print("correlation =")
print(correlation)
print("cor(x,y) =")
print(cor(x,y))
err <- abs(correlation-cor(x,y))
print("error =")
print(err)
# [1] "correlation ="
# [1] 0.0008143925
# [1] "cor(x,y) ="
# [1] 0.0008143925
# [1] "error ="
# [1] 1.12757e-17

py9. Correlation in Python

Two dice throws are independent and uncorrelated, and the calculated correlation should be near 0.

x and y correspond to 50 thousand tests of dice throw for first and second dice, respectively.

The pearsonr function from scipy.stats is used. The function cor could also have been defined with the def statement, rather than the lambda. However, lambda is usually used if the function is just one line. The name cor is used to be similar to the one in R.

# ex9.py
from __future__ import division, print_function
from numpy.random import randint
from numpy import sqrt
from pandas import DataFrame
from scipy.stats import pearsonr
cor = lambda x,y: pearsonr(x,y)[0]
num_samples = 50000
x = randint(1,7,num_samples)
y = randint(1,7,num_samples)
df = DataFrame({'x':x,'y':y})
df['a'] = df['x']-df['x'].mean()
df['b'] = df['y']-df['y'].mean()
df['ab'] = df['a']*df['b']
df['sqa'] = df['a']**2
df['sqb'] = df['b']**2
den = sqrt(df['sqa'].sum()*df['sqb'].sum())
correlation = df['ab'].sum()/den
print("correlation =", correlation)
print("cor(x,y) =", cor(x,y))
err = abs(correlation-cor(x,y))
print("error =", err)
#correlation = 0.000423927126959
#cor(x,y) = 0.000423927126959
#error = 0.0

R8. Rolling Two Dice in R

There are 5000 trials of rolling 2 dice. The sum of the 2 dices go from 2 to 12.

x1 and x2 are random samples of {1,2,3,4,5,6}, uniformly distributed. They are both vectors of length 5000.

We have to use the double brackets in the plots to return a vector from the data frame.

# ex8.R
num_samples <- 5000
x1 <- sample(6, num_samples, replace = T)
x2 <- sample(6, num_samples, replace = T)
df <- data.frame(x1,x2)
df['sum'] <- df['x1'] + df['x2']
x <- numeric(13-2)
for (i in 2:12) {
  x[i-1] <- sum(df['sum']==i)
}
df1 <- data.frame(x, row.names = 2:12)
df1['cum_x'] <- cumsum(df1['x'])
df1['cum_p'] <- df1['cum_x']/num_samples
x_lim <- c(1:6,5:1)
x_lim <- num_samples/36 * x_lim
df1['x_lim'] <- x_lim
df1['cum_lim'] <- cumsum(df1['x_lim'])
df1['cum_p_lim'] <- df1['cum_lim']/num_samples
print(df1)
err_df <- df1['cum_p_lim']-df1['cum_p']
err <- sd(err_df[,1])
print(sprintf('error = %.8f',err))
plot(2:12, df1[['cum_p']], type = 'l', col = 'red', 
     ylab = 'Cum Prob', xlab = 'x')
lines(2:12, df1[['cum_p_lim']], col = 'blue')
#     x cum_x  cum_p    x_lim   cum_lim  cum_p_lim
# 2  130   130 0.0260 138.8889  138.8889 0.02777778
# 3  305   435 0.0870 277.7778  416.6667 0.08333333
# 4  408   843 0.1686 416.6667  833.3333 0.16666667
# 5  584  1427 0.2854 555.5556 1388.8889 0.27777778
# 6  666  2093 0.4186 694.4444 2083.3333 0.41666667
# 7  851  2944 0.5888 833.3333 2916.6667 0.58333333
# 8  701  3645 0.7290 694.4444 3611.1111 0.72222222
# 9  553  4198 0.8396 555.5556 4166.6667 0.83333333
# 10 377  4575 0.9150 416.6667 4583.3333 0.91666667
# 11 301  4876 0.9752 277.7778 4861.1111 0.97222222
# 12 124  5000 1.0000 138.8889 5000.0000 1.00000000
# [1] "error = 0.00329718"

Output:

py8. Rolling Two Dice in Python

There are 5000 trials, each of rolling 2 dices. The sum of the dices is between 2 and 12.

x1 is dice1 of the 5000 rolls. x2 is dice2 of the 5000 rolls. In randint call, we give the upper number as 7 (not included). You should always look in the documentation to see if the endpoints are included which would be randint? in IPython or help(randint) in Python. For scientific work, IPython is the more productive platform.

The limiting probabilities are given in x_lim. There is 1 combination of 2 (both dices 1), 2 combination of 3 {1,2 and 2,1}, and so on. If we go above 5000 for num_samples, the error should decrease on average.

# ex8.py
from __future__ import division, print_function
import numpy as np
from numpy.random import randint
from pandas import DataFrame
import matplotlib.pyplot as plt
num_samples = 5000
x1 = randint(1,7,num_samples)
x2 = randint(1,7,num_samples)
df = DataFrame({'x1':x1, 'x2':x2})
df['sum'] = df['x1'] + df['x2']
x = np.zeros(13-2)
for i in range(2,13):
    x[i-2] = sum(df['sum']==i)
df1 = DataFrame({'x':x}, index = range(2,13))
df1['cum_x'] = df1['x'].cumsum()
df1['cum_p'] = df1['cum_x']/num_samples
x_lim = np.array([1,2,3,4,5,6,5,4,3,2,1])
x_lim = num_samples/36 * x_lim
df1['x_lim'] = x_lim
df1['cum_lim'] = df1['x_lim'].cumsum()
df1['cum_p_lim'] = df1['cum_lim']/num_samples
print(df1)
err = (df1['cum_p_lim']-df1['cum_p']).std()
print('\nerror =', err)
plt.plot(range(2,13), df1['cum_p'], 'r',
         range(2,13), df1['cum_p_lim'], 'b')
#      x  cum_x   cum_p       x_lim      cum_lim  cum_p_lim
#2   145    145  0.0290  138.888889   138.888889   0.027778
#3   272    417  0.0834  277.777778   416.666667   0.083333
#4   413    830  0.1660  416.666667   833.333333   0.166667
#5   548   1378  0.2756  555.555556  1388.888889   0.277778
#6   689   2067  0.4134  694.444444  2083.333333   0.416667
#7   837   2904  0.5808  833.333333  2916.666667   0.583333
#8   661   3565  0.7130  694.444444  3611.111111   0.722222
#9   559   4124  0.8248  555.555556  4166.666667   0.833333
#10  429   4553  0.9106  416.666667  4583.333333   0.916667
#11  303   4856  0.9712  277.777778  4861.111111   0.972222
#12  144   5000  1.0000  138.888889  5000.000000   1.000000
#
#error = 0.00353882300178

Output:

R7. Adding Mean Column and Row to a data.frame in R

The rowMeans and colMeans functions are used.

If we did not explicitly change the row.names argument in the data.frame construct, the row numbers would start with 1. This is contrast to languages such as Python, and other languages, which start at zero index.

Whenever we have new row or column indices in an assignment, a new row or column is created.

# ex7.R
x1 <- sample(0:9)
x2 <- sample(0:9)
x3 <- sample(0:9)
df <- data.frame(x1,x2,x3,row.names = c(0:9))
print('data.frame df =')
print(df)
print('First Column =')
print(df['x1'])
print('First Row =')
print(df[1,])
print('First Column, First Row = ')
print(df[1,'x1'])
print('Adding ave column')
df['ave'] <- rowMeans(df)
print('df =')
print(df)
print('Adding ave row')
df["ave_row",] <- colMeans(df)
print('df =')
print(df)
# [1] "data.frame df ="
# x1 x2 x3
# 0  2  3  4
# 1  7  8  2
# 2  1  4  5
# 3  0  1  7
# 4  9  0  1
# 5  8  6  8
# 6  3  5  6
# 7  4  7  0
# 8  6  9  3
# 9  5  2  9
# [1] "First Column ="
# x1
# 0  2
# 1  7
# 2  1
# 3  0
# 4  9
# 5  8
# 6  3
# 7  4
# 8  6
# 9  5
# [1] "First Row ="
# x1 x2 x3
# 0  2  3  4
# [1] "First Column, First Row = "
# [1] 2
# [1] "Adding ave column"
# [1] "df ="
# x1 x2 x3      ave
# 0  2  3  4 3.000000
# 1  7  8  2 5.666667
# 2  1  4  5 3.333333
# 3  0  1  7 2.666667
# 4  9  0  1 3.333333
# 5  8  6  8 7.333333
# 6  3  5  6 4.666667
# 7  4  7  0 3.666667
# 8  6  9  3 6.000000
# 9  5  2  9 5.333333
# [1] "Adding ave row"
# [1] "df ="
# x1  x2  x3      ave
# 0       2.0 3.0 4.0 3.000000
# 1       7.0 8.0 2.0 5.666667
# 2       1.0 4.0 5.0 3.333333
# 3       0.0 1.0 7.0 2.666667
# 4       9.0 0.0 1.0 3.333333
# 5       8.0 6.0 8.0 7.333333
# 6       3.0 5.0 6.0 4.666667
# 7       4.0 7.0 0.0 3.666667
# 8       6.0 9.0 3.0 6.000000
# 9       5.0 2.0 9.0 5.333333
# ave_row 4.5 4.5 4.5 4.500000

py7. Adding Mean Column and Row to a DataFrame in Python

DataFrames in pandas are similar to data.frame in R and can hold tabular data.

Every column is a specific type, such as integer or float. This is similar to R in which each column is a class such as integer or numeric.

The means function is a method of the DataFrame and the period is important here. In R, the periods are mostly treated as parts of name.

A dictionary d is created, of the 3 samples, and next a DataFrame is created based on the dictionary.

# ex7.py
from __future__ import division, print_function
from pandas import DataFrame
from numpy.random import choice
x1 = choice(10, size = 10, replace = False)
x2 = choice(10, size = 10, replace = False)
x3 = choice(10, size = 10, replace = False)
d = {'x1':x1,'x2':x2,'x3':x3}
df = DataFrame(d)
print('\nDataFrame df =\n%s' % df)
print('\nFirst Column =\n%s' % df['x1'])
print('\nFirst Row =\n%s' % df.ix[0])
print('\nFirst Column, First Row = %s' % df.loc[0,'x1'])
print('\nAdding ave column')
df['ave'] = df.mean(axis=1)
print('\ndf =\n%s' % df)
print('\nAdding ave row')
df.ix['ave_row'] = df.mean(axis=0)
print('\ndf =\n%s' % df)
#DataFrame df =
#   x1  x2  x3
#0   6   9   8
#1   3   1   6
#2   9   5   5
#3   0   3   3
#4   5   2   9
#5   2   8   7
#6   7   7   4
#7   4   0   1
#8   1   6   0
#9   8   4   2
#
#First Column =
#0    6
#1    3
#2    9
#3    0
#4    5
#5    2
#6    7
#7    4
#8    1
#9    8
#Name: x1, dtype: int32
#
#First Row =
#x1    6
#x2    9
#x3    8
#Name: 0, dtype: int32
#
#First Column, First Row = 6
#
#Adding ave column
#
#df =
#   x1  x2  x3       ave
#0   6   9   8  7.666667
#1   3   1   6  3.333333
#2   9   5   5  6.333333
#3   0   3   3  2.000000
#4   5   2   9  5.333333
#5   2   8   7  5.666667
#6   7   7   4  6.000000
#7   4   0   1  1.666667
#8   1   6   0  2.333333
#9   8   4   2  4.666667
#
#Adding ave row
#
#df =
#          x1   x2   x3       ave
#0        6.0  9.0  8.0  7.666667
#1        3.0  1.0  6.0  3.333333
#2        9.0  5.0  5.0  6.333333
#3        0.0  3.0  3.0  2.000000
#4        5.0  2.0  9.0  5.333333
#5        2.0  8.0  7.0  5.666667
#6        7.0  7.0  4.0  6.000000
#7        4.0  0.0  1.0  1.666667
#8        1.0  6.0  0.0  2.333333
#9        8.0  4.0  2.0  4.666667
#ave_row  4.5  4.5  4.5  4.500000

R6. The function sample and quantile in R

The function sample selects values randomly from a vector, with or without replacement.

With replacement, we can get a vector of length greater than the input vector. We use the unique function to see how many values are unique.

The quantile function orders our vector and gives the number corresponding to a percentage.

# ex6.R
x <- sample(0:99,size = 300, replace = TRUE)
y <- length(x)-length(unique(x))
print(sprintf('The number of repeats are %d', y))
print(sprintf('The mean is %.2f and standard deviation is %.2f',
              mean(x), sd(x)))
print('The quantiles are:')
v <- seq(.25,.75,.1)
print(quantile(x, probs = v))
## [1] "The number of repeats are 205"
## [1] "The mean is 48.77 and standard deviation is 28.79"
## [1] "The quantiles are:"
## 25%   35%   45%   55%   65%   75% 
## 24.00 34.00 42.00 54.45 65.00 73.00