x = 10; % assign x the value 10, silencing output print with ;
disp(x) % display x
x = 10 # assign x the value 10
print(x) # print the value of x
y = [1,3,5,7,9] % lists are denoted by square brackets
y(3) % calling 3rd element (MATLAB indices start at 1)
size(y) % shows that y is 1 x 5 (a row vector, by default)
length(y) % as expected, y has length 5
v = nan(2,3) % fill a 2 x 3 matrix with NaN values
size(v) % as expected, v is size (2,3)
w = repmat([1,2,3]', 2, 3) % repeats matrix twice by rows, thrice by columns
s = 1:10 % s is a list of integers from 1 to 10 inclusive
y = [1,3,5,7,9] # lists in square brackets are stored as arrays
print(y)
print(y[2]) # 3rd element (Python indices start at 0)
print(len(y)) # as expected, y has length 5
import numpy as np # NumPy: Numeric Python package
v = np.full([2,3], np.nan) # create a 2x3 matrix with NaN values
print(v)
print(v.shape) # as expected, v is size (2,3)
w = np.tile([1,2,3], (3,2)) # repeats thrice by rows, twice by columns
print(w)
s = range(10) # an iterator from 0 to 9
print([x for x in s]) # return elements using list comprehension
y = [3.14; 15; 9.26; 5]; % List with semicolons is a column vector
sum(y) % sum of all elements of y
prod(y) % product of all elements of y
max(y) % maximum value of y
min(y) % minimum value of y
range(y) % min, max value of y
mean(y) % arithmetic mean
median(y) % median
var(y) % variance
cov(y) % covar matrix = variance for single vector
corrcoef(y) % corr matrix = [1] for single vector
sort(y) % sorting in ascending order
log(y) % natural log
import numpy as np
y = [3.14, 15, 9.26, 5]
print(sum(y)) # sum of all elements of y
print(max(y)) # maximum value of y
print(min(y)) # minimum value of y
print(np.mean(y)) # arithmetic mean
print(np.median(y)) # median
print(np.var(y)) # population variance
print(np.cov(y)) # covar matrix = sample variance for single vector
print(np.corrcoef(y)) # corr matrix = [1] for single vector
print(np.sort(y)) # sort in ascending order
print(np.log(y)) # natural log
mean(y) % mean
var(y) % variance
std(y) % unbiased standard deviation, by default
skewness(y) % skewness
kurtosis(y) % kurtosis
from scipy import stats
print(np.mean(y)) # mean
print(np.var(y)) # variance
print(np.std(y, ddof = 1)) # ddof = 1 for unbiased standard deviation
print(stats.skew(y)) # skewness
print(stats.kurtosis(y, fisher = False)) # fisher = False gives Pearson definition
z = [1, 2; 3, 4] % z is a 2 x 2 matrix (Note the use of ; as row separator)
x = [1, 2] % x is a 1 x 2 matrix
%% Note: z * x is undefined since the two matrices are not conformable
z * x' % this evaluates to a 2 x 1 matrix
vertcat(z,x) % "stacking" z and x vertically
horzcat(z,x') % "stacking z and x' horizontally
%% Note: dimensions must match along the combining axis)
z = np.matrix([[1, 2], [3, 4]]) # z is a 2 x 2 matrix
x = np.matrix([1, 2]) # x is a 1 x 2 matrix
## Note: z * x is undefined since the two matrices are not conformable
print(z * np.transpose(x)) # this evaluates to a 2 x 1 matrix
b = np.concatenate((z,x), axis = 0) # "stacking" z and x vertically
print(b)
c = np.concatenate((z,np.transpose(x)), axis = 1) # "stacking" z and x horizontally
print(c)
## note: dimensions must match along the combining axis
q = -3:1:3 % specify a set of values
p = 0.1:0.1:0.9 % specify a set of probabilities
norminv(p, 0, 1) % element-wise inverse Normal quantile
tcdf(q, 4) % element-wise cdf under Student-t(4)
chi2pdf(q, 2) % element-wise pdf under Chisq(2)
%% One can also obtain pseudorandom samples from distributions
x = trnd(5, 100, 1); % Sampling 100 times from t dist with 5 df
y = normrnd(0, 1, 100, 1); % Sampling 50 times from a standard normal
%% Given sample data, we can also obtain MLE estimates of distribution parameters:
res = fitdist(x, "Normal") % Fitting x to normal dist
q = np.arange(-3,4,1) # specify a set of values, syntax arange(min, exclusive-max, step)
p = np.arange(0.1,1.0,0.1) # specify a set of probabilities
print(stats.norm.ppf(p)) # element-wise inverse Normal quantile
print(stats.t.cdf(q,4)) # element-wise cdf under Student-t(4)
print(stats.chi2.pdf(q,2)) # element-wise pdf under Chisq(2)
## One can also obtain pseudorandom samples from distributions using numpy.random
x = np.random.standard_t(df=5, size=100) # Sampling 100 times from TDist with 5 df
y = np.random.normal(size=50) # Sampling 50 times from a standard normal
## Given data, we obtain MLE estimates of parameters with stats:
res = stats.norm.fit(x) # Fitting x to normal dist
print(res) # First element is mean, second sd
x = trnd(5, 500, 1); % Create hypothetical dataset x
[h1, p1, jbstat] = jbtest(x) % Jarque-Bera test for normality
[h2, p2, lbstat] = lbqtest(x,'lags',20) % Ljung-Box test for serial correlation - Needs Econometrics Toolbox
from statsmodels.stats.diagnostic import acorr_ljungbox
x = np.random.standard_t(df=5, size=500) # Create dataset x
print(stats.jarque_bera(x)) # Jarque-Bera test - prints statistic and p-value
print(acorr_ljungbox(x, lags=20)) # Ljung-Box test - prints array of statistics and p-values
x = trnd(5, 60, 1); % Create hypothetical dataset x
subplot(1,2,1)
autocorr(x, 20) % autocorrelation for lags 1:20
subplot(1,2,2)
parcorr(x,20) % partial autocorrelation for lags 1:20
import statsmodels.api as sm
import matplotlib.pyplot as plt
y = np.random.standard_t(df = 5, size = 60) # Create hypothetical dataset y
q1 = sm.tsa.stattools.acf(y, nlags=20) # autocorrelation for lags 1:20
plt.bar(x = np.arange(1,len(q1)), height = q1[1:])
plt.show()
plt.close()
q2 = sm.tsa.stattools.pacf(y, nlags=20) # partial autocorr for lags 1:20
plt.bar(x = np.arange(1,len(q2)), height = q2[1:])
plt.show()
plt.close()
%% For loops
for i = 3:7 % iterates through [3,4,5,6,7]
i^2
end
%% If-else loops
X = 10;
if (rem(X,3) == 0)
disp("X is a multiple of 3")
else
disp("X is not a multiple of 3")
end
%% Functions (example: a simple excess kurtosis function)
%% NOTE: in MATLAB, functions can be defined in 2 locations:
%% 1) in a separate file (e.g. excess_kurtosis.m in this case) in the workspace
%% 2) in the same file as the rest of the code, BUT at the end of the file
%% function k = excess_kurtosis(x, excess)
%% if nargin == 1 % if there is only 1 argument
%% excess = 3; % set excess = 3
%% end % this is how optional param excess is set
%% m4 = mean((x-mean(x)).^4);
%% k = m4/(std(x)^4) - excess;
%% end
## For loops
for i in range(3,8): # NOTE: range(start, end), end excluded
print(i**2) # range(3,8) iterates through [3,4,5,6,7)
## If-else loops
X = 10
if X % 3 == 0:
print("X is a multiple of 3")
else:
print("X is not a multiple of 3")
## Functions (example: a simple excess kurtosis function)
def excess_kurtosis(x, excess = 3): # note: excess optional, default = 3
m4=np.mean((x-np.mean(x))**4) # note: exponentiation in Python uses **
excess_kurt=m4/(np.std(x)**4)-excess
return excess_kurt
x = np.random.standard_t(df=5,size=60) # Create hypothetical dataset x
print(excess_kurtosis(x))
y = normrnd(0, 1, 50, 1);
z = trnd(4, 50, 1);
subplot(2,2,1)
bar(y) % bar plot
title("Bar plot")
subplot(2,2,2)
plot(y) % line plot
title("Line plot")
subplot(2,2,3)
histogram(y) % histogram
title("Histogram")
subplot(2,2,4)
scatter(y,z) % scatter plot
title("Scatter plot")
y = np.random.normal(size = 50)
z = np.random.standard_t(df = 4, size = 50)
## using Matplotlib to plot bar, line, histogram and scatter plots
## subplot(a,b,c) creates a axb grid and plots the next plot in position c
plt.subplot(2,2,1)
plt.bar(range(len(y)), y);
plt.subplot(2,2,2)
plt.plot(y);
plt.subplot(2,2,3)
plt.hist(y);
plt.subplot(2,2,4)
plt.scatter(y,z);
%% Convert objects from one type to another with int8() etc
%% To check type, use isfloat(object), isinteger(object) and so on
x = 8.0;
isfloat(x)
x = int8(x);
isinteger(x)
## Convert objects from one type to another with int(), float() etc
## To check type, use type(object)
x = 8.0
print(type(x))
x = int(x)
print(type(x))