## Example Octave programs.  
##
## To run from the shell type: octave octave-examples-1.m
##
## To run from within octave, start octave then type
##    octave_examples_1 
## or 
##    source("octave_examples_1.m")
## at the command prompt.
##
## To run in Matlab, make a few modifications to the source code:
##   Change endfor, endif, endwhile, endfunction to end.
##   Change the comment character # to %.
##   Change " to '.
##   Change printf(...) to disp(sprintf(...)).
##   Change create_set to unique.
##   Change complement(A,B) to setdiff(B,A).
##   Put the function definitions in separate files.


##
## Constructing vectors.
##

## Set V to a 8-dimensional row vector with all entries equal to 0.
V = zeros(1,8);

## Set V to a 8-dimensional column vector with all entries equal to 0.
V = zeros(8,1);

## Set V to a 5-dimensional row vector with all entries equal to 1.
V = ones(1,5);

## Set V to a 5-dimensional row vector with all entries equal to -6.
V = -6 * ones(1,5);

## Set V to the 9-dimensional row vector with values 3, 4, 5, ..., 11. 
V = [3:11];

## Set V to the 9-dimensional column vector with values 3, 4, 5, ..., 11. 
V = [3:11]';

## Set V to the 5-dimensional row vector with values 3, 5, 7, 9, 11.
V = [3:2:11];

## Set V to the 3-dimensional column vector with values 3, 6, 9.
V = [3:3:11]';


##
## Arithmetic on vectors.
##

## Set V to hold the Y coordinates of the line y = 3 + 7x, where x = 1,
## 2, 3, 4, ..., 10.
V = 3 + 7 * [1:10];

## Set V to hold the Y coordinates of the quadratic function y = -1 + 2x
## - 3x^2, where x = -5, -4, ..., 4, 5.
X = [-5:5];
Y = -1 + 2*X - 3*X.^2;

## Set V to hold the Y coordinates of the logistic function exp(1-3x) /
## (1 + exp(1 - 3x)), where x = -1, -.9, ..., .9, 1.
X = [-1:.1:1];
X = exp(1 - 3*X);
V = X ./ (1 + X);

## Compute the dot product between two vectors.
V = [5:10];
W = [13:18];
D = V * W';
printf("5*13 + 6*14 + ... + 10*18 = %d\n", D);

## Another way to compute the dot product.
D = dot(V, W);
printf("5*13 + 6*14 + ... + 10*18 = %d\n", D);

## A third way to compute the dot product.
D = sum(V .* W);
printf("5*13 + 6*14 + ... + 10*18 = %d\n\n", D);


##
## Indicator vectors.
##

## Get indicator vector for the comparison V <= 4.
V = [2:8];
I = (V <= 4);
disp(V);
disp(I);

## Get the frequency of occurrence of a value greater than 2 in absolute
## value in a standard normal random vector of length 10000.
V = randn(10000,1);
p = mean(abs(V) > 2);
printf("\nFrequency(|Z| > 2) = %f\n\n", p);


##
## Subvectors and index vectors.
##

## Make a subvector of V from the index vector I that specifies
## positions 1, 3, 5, 6.  W will be the vector [2, 1, 4, 6, 6].
I = [1, 3, 6, 5, 5];
V = [2, 8, 1, 9, 6, 4, 3, 2];
W = V(I);

## Get the number of elements of V that are <= 4.
V = [2:8];
n = sum(V <= 4);
printf("\n#{V <= 4} = %d\n\n", n);

## Get the indices in V corresponding to elements that are <= 4.
V = [2:8];
I = find(V <= 4);
disp(I);
printf("\n");


##
## Set operations.
##

## U will contain the unique elements from V, in ascending order.
V = [3,2,2,3,5,1,3,8,3];
U = create_set(V);

## I is the set-theoretic intersection of A and B (unique elements, sorted).
## J is the set-theoretic union of A and B (unique elements, sorted).
## K is the set-theoretic complement of A in B (unique elements, sorted). 
A = [1,5,2,3,1];
B = [6,5,5,1,0,4];
I = intersection(A,B);
J = union(A,B);
K = complement(A,B);


##
## Sorting vectors.
##

## Sort a 10 dimensional uniform random vector in ascending order.
U = rand(10,1);
V = sort(U);
printf("Ascending sort:\n");
disp(V);

## Sort a 10 dimensional uniform random vector in descending order.
U = rand(10,1);
V = -sort(-U);
printf("\nDescending sort:\n");
disp(V);
printf("\n");

## Determine the permutation IX such that U(IX) is sorted ascending.
U = rand(10,1);
[V,IX] = sort(U);

## Get the ranks for U.
U = rand(10,1);
[V,IX] = sort(U);
[V,IX] = sort(IX);


##
## Loops and vectorization.
##

## Sum the first 100 terms of the harmonic series via a loop.
h = 0;
for i=1:100
  h = h + 1/i;
endfor
printf("1 + 1/2 + ... + 1/100 = %f\n", h);

## Vectorized sum of the first 100 terms of the harmonic series.
V = [1:100];
h = sum(1 ./ V);
printf("1 + 1/2 + ... + 1/100 = %f\n\n", h);

## Sum the first 100 terms of the alternating harmonic series via a loop.
h = 0;
for i=1:100
  h = h + (-1)^rem(i+1, 2) / i;
endfor
printf("1 - 1/2 + ... - 1/100 = %f\n", h);

## Vectorized sum of the first 100 terms of the alternating harmonic
## series.
V = [1:100]';
S = (-ones(100,1)) .^ rem(V+1, 2);
h = sum(S ./ V);
printf("1 - 1/2 + ... - 1/100 = %f\n", h);

## Generate a random integer between 10^7 and 2x10^7 and determine how
## many times it is divisible by 2.
R = 1e7 - 1 + ceil((1e7+1) * rand(1,1));

n = -1;
while ( abs(R - round(R)) < .2 )
  R = R/2;
  n = n+1;
endwhile

E = ["s."; "."];
printf("\n%d is divisble by 2 %d time%s\n\n", R * 2^(n+1), n, ...
       E(1+(n==1),:));


##
## Constructing matrices.
##

## Construct a 6x4 matrix of zeros.
M = zeros(6,4);

## Construct a 5x3 matrix of ones.
M = ones(5,3);

## Construct a 11x4 matrix of 5's.
M = 5 * ones(11,4);

## Construct a 5x3 matrix out of 3 5-dimensional column vectors.
V1 = [3, 1, 4, 0, 2]';
V2 = [1, 3, 2, 4, 0]';
V3 = [0, 4, 1, 2, 3]';
M = [V1, V2, V3];

## Construct a 2x3 matrix out of 2 3-dimensional row vectors.
V1 = [1, 3, 2];
V2 = [2, 1, 3];
M = [V1; V2];

## Construct a 5x5 matrix whose i,j entry is equal to i.
M = [1:5]' * ones(1,5);

## Construct a 5x5 matrix whose i,j entry is equal to j.
M = ones(5,1) * [1:5];

## Construct a 5x5 matrix whose i,j entry is equal to i*j.
M = [1:5]' * [1:5];

## Construct a 5x5 matrix whose i,j entry is equal to i+j.
M = ones(5,1) * [1:5] + [1:5]' * ones(1,5);


##
## Matrix slices (submatrices).
##

## Construct a 10x10 matrix out of uniform iid values.
M = rand(10,10);

## The upper left 2x2 corner of M.
C = M(1:2,1:2);

## The third column of M.
C = M(:,3);

## The fourth and fifth rows of M.
R = M(4:5,:);


##
## Elementary matrix operations.
##

## The sum and difference of two uniform random 3x5 matrices.
A = rand(3,5);
B = rand(3,5);
S = A + B;
D = A - B;

## Matrix/scalar operations.
A = rand(3,5);
B = 3 + A;
C = 3 * A;

## Construct the outer product of a multivariate normal vector with
## itself (a Wishart matrix).
V = randn(10,1);
W = V * V';

## Construct the (3x7) matrix product of a uniform random 3x5 matrix
## with a uniform random 5x7 matrix.
A = rand(3,5);
B = rand(5,7);
M = A * B;

## Construct the pointwise (Hadamard) product and the pointwise quotient
## of a uniform random 3x5 matrix with a uniform random 3x5 matrix.
A = rand(3,5);
B = rand(3,5);
P = A .* B;
Q = A ./ B;


##
## Matrix functions that act column-wise.
##

## The column sums of M (V is a 1x5 vector). 
M = rand(100,5);
V = sum(M);

## The column products of M (V is a 1x5 vector). 
M = rand(100,5);
V = prod(M);

## The column-wise maximums of M (V is a 1x5 vector). 
M = rand(100,5);
V = max(M);

## The column-wise minimums of M (V is a 1x5 vector). 
M = rand(100,5);
V = min(M);

## Sort each column of M ascending.
M = rand(20,5);
M = sort(M);


##
## Index vectors and indicator vectors for matrices.
##

## A 10x10 uniform random matrix.
M = rand(10,10);

## The indicator matrix for the elements that exceed 0.9.
Z = (M > 0.9);

## The number of elements in each column that exceed 0.9.
C = sum(M > 0.9);

## The number of elements in each row that exceed 0.9.
R = sum(M' > 0.9);

## The number of elements of M that exceed 0.9.
s = sum(sum(M > 0.9));

## The (i,j) index pairs for which M(i,j) > 0.9 is true.
[I,J] = find(M > 0.9);


##
## Centering and standardizing data.
##

## Generate a lot of data.
Z = randn(10000,100);

## Center a big array row-wise.
Z = Z - mean(Z')' * ones(1,size(Z,2));

## Center a big array column-wise.
Z = Z - ones(size(Z,1),1) * mean(Z);

## Standardize the rows of Z (i.e. linear transform so that each row has
## zero mean and unit variance).
Z = Z - mean(Z')' * ones(1,size(Z,2));
Z = Z ./ (std(Z')' * ones(1,size(Z,2)));

## Standardize the columns of Z (i.e. linear transform so that each
## column has zero mean and unit variance).
Z = Z - ones(size(Z,1),1) * mean(Z);
Z = Z ./ (ones(size(Z,1),1) * std(Z));


##
## Defining functions.
##


## The logit function.
function y = logit(x, a, b)

  u = exp(a + b*x);
  y = u / (1 + u);

endfunction


## One way to define a sample quantile: let n denote the number of
## observations, let p_k = k/n for k=1, ..., n, and let W_k be the order
## statistics. Linearly interpolate the points (p_k, W_k) to produce a
## continuous quantile function.  For p<1/n, as a special case, set the
## quantile to equal W_1.
function q = quant(V, p)

  V = sort(V);

  ## The closest number of the form i/n that is smaller than p.
  n = length(V)
  i = floor(n*p);

  ## Two special cases.
  if (i == 0)
    q = V(1);
    return;
  elseif (i == n)
    q = V(n);
    return;
  endif
  
  ## Linear interpolation.
  g = n * (p - i/n);
  q = (1-g)*V(i) + g*V(i+1);

endfunction


##
## Linear algebra.
##

## Solve the linear system Y = AX, where X and Y are 5-dimensional
## vectors and A is a 5x5 invertible matrix.
A = rand(5,5);
Y = rand(5,1);
X = A \ Y;

## Solve the linear system using matrix inversion.  The previous method
## for solving the linear system is superior.
X = inv(A) * Y;

## Singular value decomposition (SVD) of X (X = USV', U, V are unitary
## and S is diagonal).  This is the "thin" SVD where U has the same
## number of columns as X when X has more rows than columns.
X = rand(50,20);
[U,S,V] = svd(X,0);

## A Wishart matrix and its Cholesky decomposition (X = R'R, R upper
## triangular).
X = randn(5,5);
X = X * X';
R = chol(X);

## The determinant of a Wishart matrix.
X = rand(10,10);
X = X * X';
D = det(X);



##
## Further examples.
##

## Vectorized calculation of the binomial pmf.
function Q = binomial_pmf(p, n)
  
  ## The sample space.
  K = [0:n];
  
  ## The log scale pmf.
  Q = lgamma(n+1) - lgamma(K+1) - lgamma(n-K+1) + K*log(p) + (n-K)*log(1-p);
  
  ## The raw scale pmf.
  Q = exp(Q);

endfunction