I. Simplify, Expand, Factor, and Compare Expressions

Expand an expression by "epd;expr" or "epd(expr)", factor an expression by "fac;expr" or "fac(expr)", and convert a decimal to a fraction by "fra;decimal" or fra(decimal). To simplify an expression, just enter the expression.

For complex numbers, type “I” for the imaginary unit, and use the expression “abs(a+b*I)” to calculate its magnitude. To evaluate if two expressions are equivalent, enter “f(x) == g(x)” or “f(x) – g(x) == 0”, where “f(x)” and “g(x)” may be functions or number expressions. The result is “True” if the two expressions are equivalent, and “False” otherwise. Similarly, the expressions such as “x > y”, “x < y”, "x >= y", "x <= y", "f(x) > g(x)" and so on compares two numbers or expressions to determine if one is greater or less than the other.

II. Equations and Inequalities

At least three necessary elements are needed for solving an equation or inequality: The name "slv" for the operation, the expression of the equation or inequality, and the variable to be solved. The array of elements “slv; f(x); x” or “slv(f(x), x)” helps solving the equation f(x) = 0 for x. The arrays “slv; f(x) ≥ g(x); x”, “slv(f(x) ≤ g(x), x)”, “slv; f(x) > g(x); x”, and “slv; f(x) < g(x); x” find the intervals (or values) of x such that the corresponding inequality is satisfied.

In case an equation is expressed in form of f(x) = g(x), rewrite it as f(x) – g(x) = 0, and enter “slv; f(x) - g(x); x” for solution x.

Additional options can be added to the function “slv(f(x), x)”. For instance, “slv(f(x), x, C)” or “slv; f(x); x; C” finds the complex roots of f(x) = 0, where the keyword “C” is for complex roots.

III. Systems of Equations

The array "eqs(f(x, y), g(x, y))" or "eqs; f(x, y); g(x, y); x; y" helps solving a system of two equations f(x, y) = 0 and g(x, y) = 0 for the unknowns x and y, and displays the solutions in an alphabetical order.

By default, the operation “eqs” gives the solutions to all unknowns in the system. If particular variables to be solved are specified at the end such as "eqs(f(x, y), g(x, y), x, y)", other variables in the equations are treated as constants.

The array “eqs; f(x, y, z); g(x, y, z); h(x, y, z); x; y; z” or “eqs(f(x, y, z), g(x, y, z), h(x, y, z), x, y, z)” helps solving a system of three equations, where variables x, y, and z are optional unless there are additional variables in the system. Following the same pattern, use “eqs” for solving a system of four or more equations.

Example 1 Solve the linear system \(x + y + z = u + 1, 2x + y + z = 3u - 1, 3x - 2y + 2 = u,\) and \(x - y + 5z = 2u + 5\).

Example 2 Solve the nonlinear system of equations \(x^2-y^2=a\) and \(x^2+y^2=31\).

Example 3 Solve the nonlinear system by \(u^2-v^2=6\) and \(u^2-3v^2+4=0\).

I. Vector Algebra

Vector Expressions and Operations Spatial vector (vectors in R³) operations start with "vec; vec_expres", and "vec_expr" (vector expression) must be written as a linear combination of the standard basis vectors i, j, and k. The operators for vector addition, subtraction, dot product, and cross product are “+, -, *, ^ (or **)”, respectively. For simple vector operations such as "i + j", "j**k", "(2*i - j)*(4*j - 5*k)", or "x*i**(y*k)", just enter the expressions. For operations that involves several vectors, write the operation with distinct vector names, and then enter each vector name and its corresponding expression separately. For example, "u**v; u; 2*i - cos(x)*j; v; sin(x)*j - x*k" calculates the cross product of the vectors u and v. Keep in mind that the three letters "i, j, k" are reserved for the standard basis vectors, so do not use them for other vector names.

Note that “vec; expr” can carry out simple or mixed vector operations, for example, "vec;(2*i-3*k)**(-i+4*j-5*k)". For complicated vector operations involving several vectors, also use "vec; F(u, v, w); u; uexpr; v; vexpr; w; wexpr", where F(u, v, w) is the function of the vectors u, v, and w.

In a simple case, “vec; u; u; uexpr” simply returns the entries of the vector u, "vec; lnt(u); u; uexpr" gives the length (magnitude) of u, and "vec; u/lnt(u); u; uexpr" normalizes u to a unit vector.

Examples "vec; 2*u^v + a*w; u; x*i-y*j; v; cos(x)*i - y*j; w; log(y)*i + x*k" simplifies the operation 2u × v + aw among vector functions u, v, and w, where "u; x*i-y*j; v; cos(x)*i - y*j; w; log(y)*i + x*k" are the key-value pairs defining the vectors u = x*i-y*j, v = cos(x)*i - y*j, and w = log(y)*i + x*k.

These vector functions "lnt", prj", "cul", "dvg", "grd",'tnv', "ang", "hsd", and "csv" are used for finding vector length, projection, cull, divergence, gradient, tangent vector, angles between two vectors, performing the second derivative test, and checking whether a vector field is conservative.

II. Length, Angle, Projection, Area, Volume, and Orthogonal Decomposition

Calculate the length of a vector u by "lnt(u)", for example, "vec;lnt(3*i+4*j)" or "vec;lnt(u);u;cos(t)*i+sin(t)*j+2*k". The angle θ between two vectors u and v can be found by "vec;ang(uexpr, vexpr)" or "vec; ang(u, v); u; uexpr; v; vexpr", where 'ang' is the function for finding the angle between two vectors u and v, and the key-value pairs "u; uexpr; v; vexpr" define the vectors u and v. Note that the range of θ is 0 ≤ θ ≤ π. Verify the angle by "acos(u*v/(lnt(u)*lnt(v)))".

Replacing "ang" with "prj" for "vec;prj(uexpr, vexpr)", "vec; prj(u, v); u; uexpr; v; vexpr", which finds the projection of the vector u onto v, where “prj” is the function name, and the key-value pairs "u; uexpr; v; vexpr" define the two vectors. For projection, the order of u and v matters here, and "vec;prj(vexpr, uexpr)", “vec;prj(v, u); u; uexpr; v; vexpr” finds the projection of v onto u.

If v is a unit vector, calculate prj(u, v) = (u·v)v by "vec;(u*v)*v", where u·v is the component of u along v. In general, the vector u - prj(u, v) is orthogonal to v by "vec;(u-prj(u,v))*v", which is 0. Use dot product "a*b" (= 0) to verify if any two vectors a and b are orthogonal.

Calculate the area of the triangle by u and v by "vec;lnt(u**v)/2", the area of the parallelogram by u and v by "vec;lnt(u**v)", and the volume by u, v, and w by the triple scalar product "vec;abs(w*(u**v))".

III. Vector-Valued Function Calculus

One can apply the functions “lim”, “dif” and “itg” to a vector-valued function to determine its limit, derivative or integral component-wise by expressing the function as a linear combination of the standard basis vectors i, j, and k, or F₁i + F₂j + F₃k.

(1) Tangent and Normal Vectors

The array “vec;tnv(u,t); u; uexpr” or “vec;tnv(r(t), t)” finds the tangent vector of a curve parametrized by r(t) = x(t)i + y(t)j + z(t)k, where “tnv” is the operation name, and “t” is the parameter. Use the array “vec; tnv(r(t),t)/lnt(tnv(r(t),t))” to normalize the tangent vector and obtain the unit tangent vector T. Get the normal vector N by repeating the above array one more time, or use "vec; tnv(tnv(r(t),t)/lnt(tnv(r(t),t)),t)/lnt(tnv(tnv(r(t),t)/lnt(tnv(r(t),t)),t))".

One can also use "dif; r(t); t; n" or "dif(r(t), t, n)" to find the \(n\)th derivative of r(t) = x(t)i + y(t)j + z(t)k. The default n = 1 is optional. Both arrays "dif(r(t), t)" and “tnv(r(t), t)” yield the same tangent vector. It is more convenient for finding the second and third derivatives r’’(t) and r’’’(t) by "dif; r(t); t; 2" and "dif; r(t); t; 3". However, the function "tnv" returns a vector that can be used along with other vector operations.

Example Let r(t) = ti + cos(t)j + sin(t)k. The array “vec;tnv(t*i+cos(t)*j+sin(t)*k,t)/lnt(tnv(t*i+cos(t)*j+sin(t)*k,t))” finds the unit tangent vector T(t) for r(t). To find the unit normal vector N(t), apply the same operation to T(t) by “vec;tnv(tnv(t*i+cos(t)*j+sin(t)*k,t),t)/lnt(tnv(tnv(t*i+cos(t)*j+sin(t)*k,t),t))” or “tnv; T(t); t”. For the binormal vector B(t) of r(t), use “vec” to calculate T(t) × N(t) by “vec;t**n;t;i/2**(1/2)-sin(t)*j/2**(1/2)+cos(t)*k/2**(1/2);n;-cos(t)*j-sin(t)*k”, which returns \(⟨\frac{\sqrt{2}}{2},\frac{\sqrt{2}\sin t}{2},\frac{-\sqrt{2}\cos t}{2}⟩\).

(2) Tangential and Normal Components of Acceleration

To find the tangential and normal components of the acceleration for a particle moving along a path parametrized by r(t) = x(t)i + y(t)j + z(t)k, one needs to find r’(t) and r’’(t), and then apply to the formula a_T = a·T = \(\frac{\bf{a⋅v}}{||\bf{v}||}\), and a_N = a·N = \(\frac{||\bf{v × a}||}{||\bf{v}||}\), where the velocity is v(t) = r’(t) by “tnv; x(t)*i + y(t)*j + z(t)*k; t”, and the acceleration is a(t) = r’’(t) by “dif; x(t)*i + y(t)*j + z(t)*k; t; 2”.

Example let r(t) = ⟨2t, t², e^t⟩. Enter "vec;tnv(2*t*i+t**2*j+exp(t)*k,t)*tnv(tnv(2*t*i+t**2*j+exp(t)*k,t),t)/lnt(tnv(2*t*i+t**2*j+exp(t)*k,t))" for the tangential component, which is \(\frac{4t+e^{2t}}{\sqrt{4+4t^2+e^{2t}}}\). Enter "vec;lnt(tnv(tnv(2*t*i+t**2*j+exp(t)*k,t),t)**tnv(2*t*i+t**2*j+exp(t)*k,t))/lnt(tnv(2*t*i+t**2*j+exp(t)*k,t))" for the normal component, which is \(\frac{2\sqrt{4+e^{2t}+(t-1)^2e^{2t}}}{\sqrt{4+4t^2+e^{2t}}}\).

(3) Curvature

To find the curvature for a curve parametrized by r(t) = x(t)i + y(t)j + z(t)k, obtain r’(t) and r’’(t) by “dif” or “tnv” operation, and then apply the curvature formula \(κ(t) = \frac{||\boldsymbol{r}'(t) × \boldsymbol{r}''(t)||}{||\boldsymbol{r}'(t)||^3}\).

Example If r(t) = ⟨cos(t), sin(t), t⟩, "vec;lnt(tnv(cos(t)*i+sin(t)*j+t*k,t)^tnv(tnv(cos(t)*i+sin(t)*j+t*k,t),t))/lnt(tnv(cos(t)*i+sin(t)*j+t*k,t))^3" finds the curvature (κ = 0.5) of the helix r(t) for all t.

(4) Normal Vectors at Points on Space Surfaces

Combining the functions "grd" and "vec", one can find the partial derivatives r_u and r_v and the normal vector r_u × r_v at a point (u, v) on a space surface parametrized by r(u, v) = x(u, v)i + y(u, v)j + z(u, v)k.

For instance, let r(u, v) = ⟨u + v, 2u + 3v, u - v⟩. Find r_u = ⟨1, 2, 1⟩ by "vec;tnv((u+v)*i+(2*u+3*v)*j+(u-v)*k,u)", and r_v = ⟨1, 3, -1⟩ by “vec;tnv((u+v)*i+(2*u+3*v)*j+(u-v)*k,v)”. Thus, the normal vector is r_u × r_v = ⟨-5, 2, 1⟩ by “vec;(i+2*j+k)**(i+3*j-k)”.

The operation "vec; grd(F(x,y,z), x, y, z)" calculates the gradient vector for a scalar field F(x, y, z), and "vec; grd(F(x,y,z),x,f(u,v,w),y,g(u,v,w),z,h(u,v,w))*tnv(f(u,v,w)*i+g(u,v,w)*j+h(u,v,w)*k,u)" finds the derivative F_u (or F_v, F_w) by the chain rule.

(5) Curl, Divergence, Conservative and Laplacian Fields

The array “vec;cul(F, x, y, z)” calculates the curl of a vector field F, where “cul” is the function name, and F = f(x, y, z)i + g(x, y, z)j + h(x, y, z)k is the expression of a vector field, and x; y; z are the independent variables of F.

Replace “cul” with “dvg” to find the divergence of F, and replace “cul” with “csv” to determine if the field F is conservative.

(6) Properties of Curl, Divergence and Laplace Operators

Combine the functions “cul”, “dvg” and “grd” to find a Laplacian field and verify the properties of curl and divergence: (1) curl(∇f) = ∇ × (∇f) = 0; (2) dvg(curl(F)) = ∇ · (∇ × F) = 0.

I. Vectors and Matrices Manipulations

Here all matrix operations must start with the module name "mat" to distinguish it from other modules (calculus) and operations. A typical pattern for matrix or vector operation is "mat; f(A); A; matAarray", where "f(A)" represents the operation for matrix A, and the subsequent pair "A; matArray" defines the m × n matrix A with entries included in the square bracket [(a₁₁, ···, a_1n), ···, (a_m1, ···, a_mn)].

If f(A) = A, the pattern "mat; A; A; matAarray" just simply displays the entries of A. Other patterns like "mat; elf(A); A; matAarray" reduces A to an echelon form, and "mat; A*tsp(A); A; matArray" calculates the product AA^T. This basic pattern can be extended to two or more matrices for more complicated operations such as "mat; f(A, B); A; matAarray, B; matBarray", "mat; f(A, B, C); A; matAarray, B; matBarray; C; matCarray", and so on. For example, "mat; inv(P)*A*P; A; matAarray; P; matParray" gives the product P^-1AP, which are typically used for computing similar matrices or diagoalization if the columns of P have the full set of eigenvectors of A.

Row Vectors, Column Vectors, and Matrices For a row vector v = (a, b, c), "mat; v; v; [(a,b,c)]" displays v as a 1 × 3 matrix. For a column vector, "mat; v; v; [a, b, c]" displays v as a 3 × 1 matrix. The elements of a row or column are separated by a comma (,) and included in a square parenthesis. A row or column vector is a special form of a matrix. In particular, "mat; A; A; [c]" returns a 1 × 1 matrix with the only entry c.

A matrix consists of multiple rows or columns and can be entered in a similar fashion. A 2 × 2 matrix A can be entered and displayed by "mat; A; A; [(a, b), (c, d)]" with the first row "(a, b)" and the second row "(c, d)". Similarly, "mat; A; A; [(a, b, c), (x, y, z), (1, 2, 3)]" shows a 3 × 3 matrix with three rows in the order within the square bracket. The entries of each row are within a separate parathesis, and each row is seperated by a comma (,).

It is more convenient to enter a matrix row by row in order. Use the transpose function "tsp(A)" if columns are of concern.

Echelon Form and Reduced Row Echelon Form The function "elf(A)" reduces a matrix A to echelon form, and the function "rref(A)" gives the reduced row echelon form of A.

Augmented Matrices Matrices with compatible dimensions can be merged to form a new matrix. The augmented matrix [A | b] of the linear system \(A{\bf x}={\bf b}\) can be obtained by the function "mrg(A, b)" that adds \({\bf b}\) as a column to the right of A, where A is the coefficient matrix and \({\bf b}\) is the constant vector (column). With the option "1", "mrg(A,B,1)" adds the matrix B as rows under A and form a new matrix with the same number of columns.

Rotation Matrices For a given angle θ measured in radian, the functions "rtm(θ, 1)", "rtm(θ, 2)", and "rtm(θ, 3)" returns the the rotation matrices R_x, R_y, and R_y that correspond to the rotation about the x-axis, y-axis, and z-axis, respectively. They are 3 × 3 orthogonal matrices.

Jacobian Matrices In vector calculus, the Jacobian matrix is the first partial derivatives of a vector-valued function f(x, y, z) taken with respect to their parameters u, v, and w (for instance). The array "mat; jcb(a, b); a; [x(u,v,w), y(u,v,w), z(u,v,w)]; b; [u, v, w]" calculates the Jacobian matrix of f(x, y, z) = x(u, v, w)i + y(u, v, w)j + z(u, v, w)k", where "a" is the components vector "⟨x(u, v, w), y(u, v, w), z(u, v, w)⟩" of f(x, y, z), and "b" is a vector "⟨u, v, w⟩" consisting of the (symbolic) variables u, v, and w. The function "det(jcb(a,b))" calculates the determinant of the Jacobian matrix jcb(a,b).

Special Matrices The function "diag(a, b, c)" returns a diagonal matrix with diagonal entries a, b, and c, "eye(m)" returns an m × m identity matrix, "ones(m)" returns an m × m matrix with all ones, and "zeros(m) returns an m × m matrix with all zeros.

Random Matrices The function "rdm(m)" returns an m-square random matrix. In general, "rdm(m, n, type, seed, lower, upper)" returns an m × n random matrix with optional parameters "n, type, seed, lower, upper". There are 6 types of matrices indicated from 0 to 6: 0 - random matrix; 1 - diagonal matrix; 2 - upper triangular matrix; 3 - lower triangular matrix; 4 - symmetric triangular; 5 - Jordan block.

Examples "mat;rdm(4)" returns a 4 × 4 random matrix, "mat;rdm(3,2)" returns a 3 × 2 random matrix, "mat;rdm(3,3,2)" returns a 3 × 3 random upper triangular matrix, "mat;rdm(4,4,4,100,-5,5)" returns a 4 × 4 random symmetric matrix with the seed = 100 and entries randomly selected integers from -5 to 5, and "mat;rdm(4,4,5,1.5)" returns a 4 × 4 Jordan block with eigenvalue equal to 1.5.

II. Determinant, Trace, Norm, Rank, Definiteness, Diagonalizable, Nilpotent

Determinant, Trace, Norm, Rank, Condition Number For a given matrix A, the function "norm(A)" returns the Frobenius norm, "rank(A)" returns the rank, and "cnd(A)" returns the condition number.

The "norm" function has an option for different types of norms. By default, "norm(A)" returns the 2-norm of A (Frobenius). The function "norm(A,1)" returns the 1-norm of A, "norm(A,2)" returns 2-norm, and "norm(A,oo)" returns infinity norm.

For a square matrix, the function "det(A)" returns its determinant. The matrix A is singular if det(A) = 0, and nonsingular otherwise. This can be used to test if a matrix is invertiable. The "trace(A)" returns the sum of the diagonal entries of A.

True or False The function "tof(A)" (True or False) tests whether a square matrix is positive definite by default. Options can be added to the function as "tof(A, type)", and there are seven types indicated by integers from 0 to 6: 0 - diagonalizable, 1 - positive definite; 2 - negative definite; 3 - positive semidefinite; 4 - negative semidefinite; 5 - indefinite; 6 - nilpotent. To test definiteness indicated from 1 to 5, the matrix A must be symmetric or Hermitian.

Minors, Cofactors and Adjugate Matrix For the determinant (minor) of the submatrix after deleting the \(i\)th row and \(j\)th column of A, the function "mno(A, i, j)" returns the minor corresponding to the entry (i, j), and the function "cfm(A)" returns the cofactor (signed minor) matrix of A (for all entries). The transpose of the cofactor matrix of A is the ajugate (classic adjoint) of A. Note that the minor corresponds to the entry (1, 1) of A is obtained by mno(A,0,0), rather than mno(A,1,1).

III. Matrix and Vector Operations

Addition, Subtraction, Multiplication, and Power The operators +, -, * also apply to matrix addition, subtraction, scalar and matrix multiplication. For multiplying two matrices A and B, "A*B" and "B*A" calculate the products of AB and BA, respectively, provided they are defined. Enter "A**n" or "A^n" for the power of a matrix A, where "A" must be a square matrix and "n" is a whole number exponent.

For example, the "A+3*B" gets the sum of A + 3B, "A**2-2*A*B + B**2" returns the result of A² - 2AB + B².

Transpose, Conjugate, and Conjugate Transpose: The function "tsp(A)" returns A^T, the transpose of A. For complex matrices, "tsp(A)" returns the conjugate transpose of A, and "cjg(A)" just returns the conjugate of A.

Matrix Inverse Both "inv(A)" and "A^(-1)" return the inverse of a nonsingular matrix A. If A is not square, "inv(A)" returns the Moore-Penrose pseudoinverse, but "A^(-1)" is not defined.

Note that matrix power with integer exponents only applies to square matrices. The result of "A^(-2)" can be understood as the squared inverse of A. Thus, "A*B**(-1)" gives the result for AB^-1, and "3*A*B**3+4*C**(-2)" for 3AB³ + 4C^-2, provided they exist.

Matrix Polynomials A matrix polynomial is a polynomial with square matrices as variables. Given an ordinary scalar-valued polynomial \(f(t)=a_0+a_1t+a_2t^2+\cdots+a_nt^n\), the matrix polynomial is \(f(A)=a_0I+a_1A+a_2A^2+\cdots+a_nA^n\). For example, "A^2-3*A+2*eye(2)" calculates the matrix polynomial \(f(A)=A^2-3A+2I_2\).

IV. Linear System, Linear Independence, and Vector Coordinates

Solve a Linear System Many problems in linear algebra reduce to finding the solutions of a system of linear equations, and solving a linear system is an important part of linear algebra. To solve a nonhomogeneous system of linear equations Ax = b, form the augmented matrix M = [A, b], and use the function "rref(M)" to obtain the solution, where "A" represents the coefficients matrix and "b" the constant vectors (column). The same results can be obtained by "rref(mrg(A,b))", where "mrg(A,b)" adds b as a column to the right of A. If the system has a unique solution, the solution is in the rightmost column of "rref(M)" or "rref(mrg(A,b))".

If there is a "0 = 1" row in rref(mrg(A,b)), the system has no solution. If there is a free variable in A, the system has infinitely many solutions. The corresponding homogeneous linear system Ax = 0 always has a zero solution. Use "rref(A)" to solve the homogeneous system. It has the zero solution if rref(A) is an identity matrix. If there is a free variable, the system has nonzero solutions. The number of free variables is the dimension of null space of A. The general solution of Ax = b can be expressed as the sum of a particular solution of Ax = b and the null space of A.

Determine Linear Independence To show a set of vectors S = \(\{{\bf v}_1,{\bf v}_2,\cdots,{\bf v}_n\}\) is linearly independent, form a matrix A = [\({\bf v}_1,{\bf v}_2,\cdots,{\bf v}_n\)] whose columns are the given vectors in S, and reduces A to echelon form by by "elf(A)". The set is linearly dependent if there is a zero row in elf(A), and independent otherwise.

Or solve the homogeneous system Ax = 0 by "rref(A)". The set is linearly dependent if rref(A) has a free column, and independent otherwise.

Determine a Vector in a Span Given a spanning set S = \(\{{\bf v}_1,{\bf v}_2,\cdots,{\bf v}_n\}\), to determine if a vector b is in span(S), solve the linear system Ax = b by "rref(mrg(A,b))", where the columns of A are the vectors in S. The b is not in span(S) if rref(mrg(A,b)) has a row "0 = 1", which means the system does not have a solution.

Find Coordinate Vectors Given a basis S = \(\{{\bf v}_1,{\bf v}_2,\cdots,{\bf v}_n\}\) of a vector space V, to obtain the coordinate vector of a vector b ∈ V relative to S, solve the linear system Ax = b by "rref(mrg(A,b))", where the columns of A are the vectors in S. The unique solution is the coordinate vector [b]_S.

Obtain a Vector Orthogonal to a Subspace Given a set S = \(\{{\bf v}_1,{\bf v}_2,\cdots,{\bf v}_n\}\) of vectors, to obtain a vector x that is orthogonal to the vectors in S, solve the linear system Ax = 0 by "rref(A)", where the rows of A are the vectors in the set. This is because the null space of A is perpendicular to the row space of A.

Obtain Eigenvectors For a given n-square matrix A, to obtain the eigenvectors (or a basis of the eigenspace) belonging to an eigenvalue λ, solve the linear system (A - λI)x = 0 by "rref(A-λ*eye(n))", where "eye(n)" gives an n × n the identity matrix.

Least Squares Method is used for a linear system \(A{\bf x}={\bf b}\) that has more equations than unknowns and the system has no solutions. The method is used for solving \(A^TA\hat{\bf x}=A^T{\bf b}\) for the best solution \(\hat{\bf x}\). If \(A\) has \(n\) independent columns, \(A^TA\) is invertible, and \(\hat{\bf x}=(A^TA)^{-1}A^T{\bf b}\) by "inv(tsp(A)*A)*tsp(A)*b" minimizes the errors \(||{\bf b}-A{\bf x}||\), where \(A\hat{\bf x}=A(A^TA)^{-1}A{\bf b}\) is the projection of \({\bf b}\) onto the column space of \(A\).

V. Vector Space, Dimension, and Basis of Subspaces

Four Fundamental Subspaces For an m × n matrix A, the function "nsp(A)" returns a list of row vectors (or basis) that span the null space of A, and "nsp(tsp(A))" returns a list of vectors (row) that span the null space of A^T (left null space), where "nsp" (null space) is the function for finding a basis of the null space (or kernel of a linear mapping with matrix A).

The function "rsp(A)" returns a list of vectors (row) that span the row space of A, and "rsp(tsp(A))" returns a list of vectors (row) that span the column space of A, where "rsp" (row space) is the function for finding a basis for a row space of A.

Also use "rref(A)" to obtain the dimension and a basis for each subspace. The nonzero rows of rref(A) form a basis for the row space, and the number of the nonzero rows is the dimension of the row space or the rank of the matrix A. The columns in A that correspond to the pivot columns in rref(A) form a basis of the column space. If the number of pivots is r, the dimension of N(A) is n - r, and the dimension of N(A^T) is m - r.

Dimension and Basis of a Vector Space Given a set of vectors \(S=\{{\bf v}_1,{\bf v}_2,\cdots,{\bf v}_n\}\) that spans a subspace W, form a matrix A whose rows are the vectors in S, and evaluate the echelon matrix "elf(A)" to determine the dimension and a basis of W. The number of the nonzero rows in echelon form is the dimension of W, and the nonzero rows form a basis of W.

Kernel and Image of Linear or Matrix Mapping Since an \(m×n\) matrix \(A\) can be viewed as a linear mapping, the image of the mapping is the column space of \(A\) by "rsp(tsp(A))", and the kernel of the mapping is the null space of \(A\) by "nsp(A)".

Rank and Nullity of Linear or Matrix Mapping If A is the matrix for a linear mapping \({\bf T}:V→U\), the rank of the mapping is the dimension of the image of \({\bf T}\) (or the column space of A), and the nullity of the mapping is the dimension of kernel of \({\bf T}\) (or the null space of A). Reduce the matrix A to echelon form by "elf(A)". The number of pivots is the rank of \({\bf T}\), the number of non-pivot columns is the nullity of \({\bf T}\), and the sum of the rank and nullity is the dimension of \(V\).

VI. Linear Mappings, Matrix Representations, and Change of Basis

Determine the Matrix for a Linear Operator Given a linear operator \({\bf T}:R^n→R^n\), there is a matrix A such that \({\bf T}({\bf v})=A{\bf v}\) relative to the standard basis of \(R^n\). If \(S=\{{\bf v}_1,{\bf v}_2,\cdots,{\bf v}_n\}\) is a non-standard basis of \(R^n\), the columns of the matrix \([{\bf T}]_S\) for \({\bf T}\) relative to the basis \(S\) are the solution to each linear system \(B{\bf x}={\bf T}({\bf v}_i)\), where the columns of \(B\) are the basis vectors in \(S\).

Instead of solving each linear system one at a time, use "rref(mrg(B,A*B))" to get the matrix for \({\bf T}\), where \(A\) is the matrix for \({\bf T}\) relative to the standard basis, the columns of \(B\) are the basis vectors in S, and the matrix for \({\bf T}\) relative to the basis \(S\) is on the right of the reduced row echelon form. The main concern is to choose a basis that gives the best matrix (diagonal) for \({\bf T}\).

Matrix for General Linear Transformation Any linear mapping \({\bf T}:R^n→R^m\) is assigned an \(m×n\) matrix \(A\) by \({\bf T}({\bf v})=A{\bf v}\) relative to the standard bases of \(R^n\) and \(R^m\). But both spaces \(R^n\) and \(R^m\) have non-standard bases. If \(V=\{{\bf v}_1,{\bf v}_2,\cdots,{\bf v}_n\}\) is a basis for \(R^n\) and \(W=\{{\bf w}_1,{\bf w}_2,\cdots,{\bf w}_m\}\) is a basis for \(R^m\). The matrix for \({\bf T}\) relative to the bases \(V\) and \(W\) can be obtained by "rref(mrg(W, A*V))", where \(A\) is the matrix for \({\bf T}\) relative to the standard bases, the columns of \(V\) are the basis vectors of \(R^n\), the columns of \(W\) are the basis vectors of \(R^m\), and the matrix \([{\bf T}]_{V,W}\) for \({\bf T}\) relative to the bases \(V\) and \(W\) is on the right of the reduced row echelon form.

Change of Basis Matrix Given two different bases \(V=\{{\bf v}_1,\cdots,{\bf v}_n\}\) and \(W=\{{\bf w}_1,\cdots,{\bf w}_n\}\) of \(R^n\), the change of basis matrix from \(V\) to \(W\) is \(P=[[{\bf w}_1]_V,\cdots,[{\bf w}_n]_V]\), where the columns \([{\bf w}_i]_V\) are the coordinate vectors of \({\bf w}_i\) relative to the basis \(V\) for each \(i\). This means to express each basis vector \({\bf w}_i\) in the basis of \(V\), and the question equals to find the coordinate vector \([{\bf w}_i]_V\) by solving the linear system \(V{\bf x}={\bf w}_i\). Thus, the change of basis matrix \(P\) can be obtained by "ref(mrg(V, W))", where the columns of matrix \(V\) are the basis vectors in \(V\) and the columns of \(W\) are the basis vectors in \(W\). If \(V\) is the standard basis of \(R^n\), then \(P=[{\bf w}_1,\cdots,{\bf w}_n] = W\).

The change of basis matrix \(Q\) from \(W\) to \(V\) is \(Q=P^{-1}=[[{\bf v}_1]_W,\cdots,[{\bf v}_n]_W]\). It can be obtained by "rref(mrg(W, V))", or find the inverse of \(P\) directly by "ivs(P)".

VII. Inner Product Space, Projection, and Orthogonality

Inner Product on \(R^n\) or \(C^n\) For vectors u and v in \(C^n\), "tsp(v)*u" returns the inner product of ⟨u, v⟩, and "tsp(u)*v" returns the inner product of ⟨v, u⟩. Note that the results of "tsp(v)*u" and "tsp(u)*v" are a single-entry matrix instead of a scalar.

The function "dot(u,v)" or "dot(v,u)" returns a scalar of the dot product ⟨u, v⟩ or ⟨v, u⟩. For complex vectors, "tsp(u)" gives the conjugate transpose of u, and "dot(tsp(u),v)" or "dot(u,tsp(v))" calculates the dot product of complex vectors.

Projection of v onto w If w is nonzero, the component (Fourier coefficient) \(c\) of v along w can be calculated by "dot(v,w)/dot(w,w)", and the projection of v along w is \(c\)w by "w*dot(v,w)/dot(w,w)".

Projection of v along a Subspace W Suppose W = span\(\{{\bf w}_1,{\bf w}_2,\cdots,{\bf w}_k\}\) is a subspace of \(R^n\). The projection of a vector \({\bf v}∈R^n\) in W equals \(A(A^TA)^{-1}A^T{\bf v}\) by "A*inv(tsp(A)*A)*tsp(A)*v", where the columns of \(A\) are the the vectors \({\bf w}_1,{\bf w}_2,\cdots,{\bf w}_k\). The matrix \(A(A^TA)^{-1}A^T\) is called projection matrix. If W = \({\bf w}\) is a nonzero vector, the projection matrix is \(\frac{{\bf w}{\bf w}^T }{{\bf w}^T{\bf w} }\).

Gram-Schmidt Orthogonalization Process Given a basis \(S=\{{\bf v}_1,{\bf v}_2,\cdots,{\bf v}_n\}\), the function "grs(A)" returns a set of orthogonal vectors from \(S\), where the rows of \(A\) consists the vectors in \(S\) and the rows of grs(A) are orthogonal vectors corresponding to \(S\). The function "grs(A,1)" further normalizes the orthogonal vectors to form an orthonormal set.

Verify Orthogonal Matrix Check if the result of "inv(P)*P" is an identity matrix. If it is an identity matrix, \(P\) is orthogonal. To verify if \(P\) is orthogonal is to verify that each row or column has unit length and that each pair of rows or columns is orthogonal.

VIII. Eigenvalues/Eigenvectors, Singular Values/Vectors, Quadratic Forms, and Diagonalization

Characteristic Polynomials The function "chp(A)" returns the characteristic polynomial \(p(λ)\) of a matrix \(A\), which is defined by det\((A-λI)\).

Eigenvalues and Eigenvectors The function "eig(A)" returns the eigenvalues and eigenvectors of a matrix \(A\). The output of "eig(A)" (eigenvalues or eigenvectors) displays each eigenvalue, its algebraic multiplicity and the corresponding eigenvectors.

The eigenvalues of \(A\) are the roots of the characteristic polynomial \(f(λ)=\) det\((A-λI)\) and can be obtained by solving the equation det\((A-λI)\) = 0 using "slv;f(λ);λ". The eigenvector with eigenvalue \(λ\) can be obtained by solving the linear system \((A-λI){\bf x}={\bf 0}\) using the "rref" function.

Normal Matrices \(A\) is normal if \(A\) commutes to its transpose. To verify if a square matrix \(A\) is normal is to verify the equality \(AA^{T}=A^{T}A\) by "A*tsp(A) - tsp(A)*A". If the result is a zero matrix, then \(A\) is normal. Hermitian (symmetric), skew-Hermitian (skew-symmetric), positive (definite) matrices are normal matrices.

Diagonalization All normal matrices (or linear operators on a complex finite-dimensional inner product space) are diagonalizable. To determine if a matrix is diagonalizable, use "tof(A,0)" and it returns "True" or "False", where the option 0 in the function "tof" is used to test whether a square matrix is diagonalizable.

The function "dgl(A)" returns two matrices \(P\) and \(D\) such that \(D=P^{-1}AP\) and \(D\) is diagonal, where the function "dgl" (diagonalization) has two optional parameters dgl(A, normalize, PD). By default, "dlg(A)" returns two matrices \(P,D\) and the columns of \(P\) are not normalized (normalize = 0). If normalize = 1, the columns of \(P\) are normalized. The function returns the matrix \(P\) if PD = 1, and it return \(D\) if PD = 2. For example, the array "mat; inv(dgl(A,0,1))*A*dgl(A,0,1); A; matAarray" verifies the relation \(P^{-1}AP=D\).

The diagonal matrix \(D\) is similar to \(A\). If the n-square matrix \(A\) has a full set of eigenvectors, choose the eigenvectors as the basis of \(R^n\) by "eig(A)", form a matrix \(P\) with the eigenvectors as its columns, and then the diagonal entries of \(D=P^{-1}AP\) are the eigenvalues of \(A\).

Diagonalizing Quadratic Forms A quadratic form determines a symmetric matrix \(A\). Use the function "dgl(A)" to find the orthogonal matrix \(P\) such that \(D=P^TAP\). Check the example "mat;dgl(A);A;[(8,3),(3,2)]".

Matrix Exponential and Power If \(A=PDP^{-1},e^{At}=I+A+\frac{(At)^2}{2!}+\cdots+\frac{(At)^n}{n!}+\cdots=Pe^{Dt}P^{-1}\), where the diagonals of \(D\) are the eigenvalues of \(A\), and the columns of \(P\) are the full set of eigenvectors of \(A\). For example, A = [(1,2),(1,0)], expt(A) can be obtained by "eig(A)", and "mat;P*D*inv(P);P;[(2,1),(-1,1)];D;exp(diag(2,-1))" for t = 1. The matrix power \(A^k=XD^kA^{-1}\) can be obtained by "PD^kP^(-1)".

Singular Values and Singular Vectors Any \(m×n\) matrix \(A\) can be factorized as \(A=UΣV^T\) (singular value decomposition), where \(U=[{\bf u}_1,\cdots,{\bf u}_r], V=[{\bf v}_1,\cdots,{\bf v}_r]\) are orthogonal matrices consisting of orthonormal vectors (called singular vectors), and \(Σ=\) diag\((σ_1,\cdots,σ_r)\) for singular values \(σ_1≥σ_2≥\cdots≥σ_r>0\). The singular values \(σ_i\) and singular vectors \({\bf v}_i\) can be obtained by "eig(tsp(A)*A)" because \(A^TA{\bf v}_i=σ_i^2{\bf v}_i\). Since \(A[{\bf v}_1,\cdots,{\bf v}_r]=[{\bf u}_1,\cdots,{\bf u}_r]Σ\) or \(A{\bf v}_i=σ_i{\bf u}_i\), the singular vectors \({\bf u}_i=\frac{A{\bf v}_i}{σ_i}\).

Difference Equations In a difference equation \({\bf u}_{k+1}=A{\bf u}_k\), each step multiplies by \(A\), and the solution is \({\bf u}_k=A^k{\bf u}_0=PD^kP^{-1}{\bf u}_0\), where the columns of \(P\) are the full set of eigenvectors \(\{{\bf x}_1,\cdots,{\bf x}_n\}\) of \(A\). If \({\bf u}_0=c_1{\bf x}_1+\cdots+c_n{\bf x}_n,{\bf u}_k=c_1λ_1^k{\bf x}_1+\cdots+c_nλ_n^k{\bf x}_n\). For example, the Fibonacci sequence \({\bf u}_k=[F_{k+1},F_k]^T=A^k{\bf u}_0\) can be obtained by "mat;a*b*inv(a)*u;a;[((1-sqrt(5))/2,(1+sqrt(5))/2),(1,1)];b;diag((1-sqrt(5))/2,(1+sqrt(5))/2);u;[0,1]", where \(A\) = [(1,1),(1,0)], and the eigenvalue and eigenvectors can be obtained by "mat;eig(A);A;[(1,1),(1,0)]".

System of Differential Equations Extend a 1 by 1 single differential equation (DE) \(u'=u\) or \(u'=λu\) to a system of n by n linear constant-coefficient DEs \(\frac{d{\bf u}}{dt}=A{\bf u}\) or \({\bf u}'=A{\bf u}\). The general solution is \({\bf u}(t)=c_1e^{λ_1t}{\bf x}_1+\cdots+c_ne^{λ_nt}{\bf x}_n\), where \(λ_1,\cdots,λ_n\) and \(\{{\bf x}_1,\cdots,{\bf x}_n\}\) are the eigenvalues and the full set of eigenvectors of \(A\) by "eig(A)". For example, the system \(x'=y+z,y'=x+z,z'=x+y\) with \({\bf u}=[x,y,z]^T, A\) = [(0,1,1),(1,0,1),(1,1,0)]. The general solution is \({\bf u}(t)=c_1e^{-1}[-1,1,0]^T+c_2e^t[-1,0,1]^T+c_3e^{2t}[1,1,1]^T\). If \({\bf u}_0=[2,3,1]^T, [c_1,c_2,c_3]=[1,-1,2]\) by "mat;rref(mrg(P,u));P;[(-1,-1,1),(1,0,1),(0,1,1)];u;[2,3,1]".

IX. Matrix Decomposition

Common matrix decomposition methods are implemented. They are (1) QR-decomposition ("qrd"), (2) LU-decomposition ("lud"), (3) singular value decomposition ("svd"), (4) Cholesky decomposition ("chl"), (5) rank decomposition ("rkd"), and (6) LDL-decomposition ("ldl").

The function "qrd(A)" factors A as A = QR; "lud(A)" factors A as A = LU. In a similar fashion, the function "ldl(A)" returns matrices L and D for which A = LDL^*; "cho(A)" returns matrix L for A = LL\(^*\); "svd(A)" returns matrices U, Σ and V for A = UΣV; "rkd(A)" returns matrices C and F for A = CF.

Modules starting with "plt" are for data visualization and exploration (e.g., histograms, scatter plots, boxplots, and graphs of cdfs, pdfs, and ecdfs); modules starting with "stt" are for symbolic computing, data analysis and hypothesis testing (e.g., mean, variance, probabilities, statistical models, hypothesis testing, and random processes).

I. Common Distributions

1. Discrete distributions Discrete distribution modules include Bernoulli (ber), binomial (bin), Poisson (poi), geometric (geo), negative binomial (nbn), hypergeometric (hyp), and discrete uniform (die). Check stt;ber(0.6);v || stt;poi(10);std | stt;bin(10,0.5);P(x>2) || stt;geo(0.7);m || stt;nbn(8,0.6);m || stt;hyp(20,5,12);rvs || stt;die(10);cdf(x) ||. Note that you need to provide appropriate parameter(s) for each discrete distribution. Some distributions have one parameter, and some have two or more parameters.

2. Continuous distributions Continuous distribution modules include normal (n), exponential (epn), uniform (uni), Gamma (gam), Beta (bet), Cauchy (cau), F-distribution (f), Student-t distribution (t), Chi-square (chi). Check stt;n(3,1);q(0.8) || stt;epn(0.5);svf(1.2) || stt;uni(0,2);V(x) || stt;gam(2,3);P(x>2&x<4) || stt;bet(0.5,2);E(x) || stt;cau(0,1);cdf(x) || stt;f(8,11);P(x>13) || stt;t(2);P(x>2@x<-2) || stt;chi(1);q(0.975) ||.

3. User defined distributions Use "drv" (discrete random variable) to create a discrete random variable (r.v.). For example, stt;drv(0.75^x/3);P(x>2) creates a geometric random variable with probability mass function (pdf) \(f(x)=\frac{1}{3}(\frac{3}{4})^x\) for \(x=1,2,\cdots\), and \(\sum_{x=1}^∞f(x)=1\). For a finite discrete random variable, list the probability mass function in two different rows. For example, stt;drv((0,1,2),(2/4,1/4,1/4));cdf(x) ||. The first row lists the possible values the random variable can take, and the second lists their corresponding probabilities that must sum to 1.

Use "crv" (continuous random variable) to create a continuous random variable with specified pdf and interval (range). In this example "stt;crv(3*x^2/26,1,3);E(x)", the pdf is \(f(x)=\frac{3x^2}{26}\) and the range or interval the r.v. can take is [1, 3]. Note that the pdf must be integrated to 1 over this interval.

4. Compute probabilities Use the function "P(.)" for computing probabilities. For example, "P(2)" standards for the probability \(P(X=2)\), "P(x!=2)" for \(P(X≠2)\), "P(x>2)" for \(P(X>2)\), and "P(x<=2)" for \(P(X≤2)\). For compound inequalities, "P(y>2&y<=5)" or "P(2< y<=5)" stands for the probability \(P(2< Y≤5)\) or \(P(Y>2\text{ and }Y≤5)\), and "P(z>=6@z<3)" standards for the probability \(P(Z≥6∪Z<3)\) or \(P(Z≥6\text{ or }Z<3)\). Use similar combinations to compute probabilities. Check stt;bin(10,0.5);P(x>=6) || stt;epn(0.8);P(x<1) || stt;geo(0.7);P(y!=2) || stt;n(0,1);P(x>1.96@x<-1.96) || stt;n(0,1);P(x>-1&x<1) || stt;poi(6);P(2@5) || stt;bin(10,0.5);P(3<=w<=7) ||.

Use the function "P(.|.)" to compute conditional probabilities. For example, "P(x>5|x>=2)" represents the conditional probability \(P(X>5|X≥2)\). Similarly, "P(x<=4|x<8)" represents the probability \(P(X≤4|X<8)\). Check stt;poi(8);P(x<=4|x<8) || stt;poi(8);P(x>4|x>2) || stt;epn(1.2);P(x>2|x>1) || stt;geo(0.6);P(x>5|x>4) ||.

5. Compute mean, variance, skewness, kurtosis, entropy and moments Use the letter "m, v, s, k, h" to obtain the mean, variance, skewness, kurtosis, and entropy of the r.v., respectively. Use "std, med" to get the standard deviation and median of the r.v., respectively. Use "M(r)" for computing moments, where "r" is a positive integer representing the r-th moment. Check stt;ber(0.6);h || stt;bin(20,0.4);m || stt;poi(9);med || stt;n(3,1.5);s || stt;epn(2);k || stt;gam(4,3);std || stt;uni(1,4);v || stt;poi(3);M(3) || stt;n(4,2);M(4) || stt;epn(4);s ||.

You can obtain the above quantities by "E(x), V(x), s(x), k(x), H(x), std(x), med(x)", which calculate mean, variance, skewness, kurtosis, entropy, standard deviation and median, respectively. You can also compute the means and variances of linear combinations of the r.v.. For example, "E(2*x+1)", "V(2x+1)" compute the mean \(E(2X+1)\) and and \(\text{Var}(2X+1)\). In particular, "E(exp(t*x));x" returns the moment generating function for the r.v. ("x"), "E(exp(w*I*x));x" give the characteristic function of the r.v. ("x"). Ensure to add "x" to the end to tell the program that letter "x" represents the r.v. (instead of others letters (like "t", "w" or "I") in the expression). For example, stt;n(0,1);E(2*x+1) || stt;n(2,1.5);V(2*x+1) || stt;n(0,1);E(exp(t*x));x || stt;epn(3);E(exp(t*x));x || stt;n(0,1);E(exp(w*x*I));x || stt;crv(3*x^2/26,1,3);E(x) || stt;crv(3*x^2/26,1,3);std(x) ||.

6. Compute cdf, pdf, percentile, and (inverse) survival function Use "pdf(a), cdf(a), svf(a), isf(p), q(p)" to compute the value of pdf, cdf, survival function, inverse survival function, and the p-th percentile (quantile), respectively, for a given number "a" or "p" (for p ∈ [0, 1]). For example, stt;bin(10,0.7);cdf(6) || stt;poi(12);pdf(10) || stt;t(8);q(0.95) || stt;epn(0.8);svf(1) || stt;epn(2);isf(0.6) ||.

Use "pdf(x)" and "cdf(x)" to show the formulas for the pdf and cdf of the r.v. "x". For example, stt;n(0,1);pdf(x) || stt;epn(1);cdf(x) ||.

7. Simulating samples Use "rvs(n)" for simulating a random sample of size "n" from common distributions or user-defined distributions. For example, stt;chi(1);rvs(3) || stt;t(8);rvs(20) || stt;hyp(15,4,8);rvs(20) ||.

Use keyword "d" to obtain some summary statistics from simulated samples. Observe how these statistics change as the sample size n increases and examine the relationship between the sample statistics and population parameters. For example, stt;poi(20);rvs(30);d || stt;poi(20);rvs(300);d ||.

II. Descriptive Statistics

For a data set with one row of numbers, you can get the following statistics by keyword "d" (descriptive): sample size, min, max, mean, standard deviation (std), variance, skewness, kurtosis, median, inter-quantile range, and the standard error for mean estimate. Check stt;dat(5,1,10,5,6,10,9,9,4,8);d || stt;dat(1.99,-0.29,-0.93,-0.82,1.73,-0.57,-0.91,0.2,-0.42,0.2);d ||.

Use keyword "z" for the z score of each value in the sample, relative to the sample mean and standard deviation. For example, stt;dat(11.73,13.93,12.79, 13.92,12.34,11.15,13.14,11.86,16.61,13.2);z ||.

Compute the p-th quantile or 100p% percentile by function "q(p)", where "p" is a number between 0 and 1. For example, stt;dat(5,1,10,5,6,10,9,9,4,8);q(0.6) ||. You can also compute the empirical cumulative density function for any number "a" by the function "ecdf(a)". For example, stt;dat(5,1,10,5,6,10,9,9,4,8);ecdf(2) || stt;dat(5,1,10,5,6,10,9,9,4,8);ecdf(6.7) ||.

III. Data Visualization and Exploration

To obtain a graph, start with "plt" rather than "stt". Use the form "plt;dat(one row);ecdf" to obtain the graph of the empirical cumulative density function, use form "plt;dat(one rwo);kde" to obtain the graph of Gaussian kernel density estimate for a continuous random variable, use "plt;dat(one rwo);pt" or "plt;dat((rowX),(rowY));pt" to obtain the scatter plot of your data (one or two rows for variables X and Y), and use form "plt;dat((row1),(row2),...,(rowk));box" to obtain the boxplot for one or multiple rows of your data. Check the examples: plt;dat((-1,2,1,0.5,4,3),(-2.8,3,2,0.8,7,5));pt || plt;dat(3,12,5,5,6,9,7,7,5,6,10,12,6,7,8,7,8,10,7,6);ecdf || plt;dat(4.97,4.96,6.38,4.85,5.89,5.13,6.38);box || plt;dat(3,12,5,5,6,9,7,7,5,6,10,12,6,7,8,7,8,10,7,6);kde || plt;dat((3,12,5,5,6,9,7,7,5,6,10,12,6),(2,8,5,9,-5,4,7,7,8,7,8,10,7));box ||.

You can plot the probability density function (pdf), cumulative density function (cdf) and histogram (hist) of samples randomly generated from some common distributions. For example, plt;n(2,1);pdf || plt;n(2,1);cdf || plt;n(2,1);hist || plt;poi(5);cdf || plt;poi(5);hist ||. Note that "n(2,1)" standards a normal distribution with mean 2 and standard deviation 1, and "poi(5)" for a Poisson distribution with mean 5.

IV. Hypothesis Testing

1. Student-t test Use "t1" for one-sample t-test, "t2" for two-samples t-test, and "pt" for paired-samples t-test. For example, stt;dat(2.84,3.76,2.61,4.29,2.15,1.98,3.21,4.38,2.9,3.11,2.89,2.55,3.74,3.93,2.95);t1;3 || tests the mean equals 3; stt;dat((4,2,1,5,0,5,1,1,3,3,0,7,2,3,7),(6,2,6,5,2,0,4,3,8,1,3,7,4,5,5));t2 || tests if the two means are significantly different; stt;dat((5,2,8,4,1),(7,3,6,5,3));pt || tests if the two means are significantly different before and after some treatment.

2. Goodness of fit test Use "chi" to determine whether there is a statistically significant difference between the expected frequencies and the observed frequencies in one or more categories of a contingency table (test of homogeneity). For example, stt;dat((5,8,9,2,1,7),(6,6,6,4,4,6));chi ||, where the first row is the observed counts and the second row gives the expected frequencies. When just observed frequencies are given, it is assumed that the expected frequencies are uniform and given by the mean of the observed frequencies. For example, stt;dat(5,8,9,2,1,7);chi ||.

Use "n" to test normality. Use "ks1" (Kolmogorov-Smirnov test one sample) to test if a sample comes from a certain population distribution, and use "ks2" (Kolmogorov-Smirnov test one sample) to test if two samples come from the same distribution. For example, stt;dat(6.03,5.41,4.01,6.17,3.93,5.24,5.97,6.4,4.58,6.11,5.43,5.65,6,4.31,4.42);n || stt;dat((1,7,9,7,12,9,12,9,6,12),(3,7,7,8,9,10,8,7,6,13));ks2 ||. To use "ks1", specify the distribution by keyword such as "n", "t", "f", "chi", "exp", "poi", "cau", "gam", "bet", and "uni", followed by its necessary parameter(s). For example, stt;dat(6.03,5.41,4.01,6.17,3.93,5.24,5.97,6.4,4.58,6.11,5.43,5.65,6,4.31,4.42);ks1;n;3;1 || stt;dat(6.03,5.41,4.01,6.17,3.93,5.24,5.97, 6.4,4.58,6.11,5.43,5.65,6,4.31,4.42);ks1;n;5;1 || stt;dat(12,4,5,9,8,10,5,9,1,5,9,9,6,7,3);ks1;poi;7 || stt;dat(1.27,2.3,1.56,0.48,4.44,1.52,1.09,0.19, 1.5,1.27,2.2,0.59,2.43,1.67,0.45);ks1;epn;0.8 ||.

3. Analysis of variance (ANOVA) Use "f1" for one-way ANOVA and "f2" for two-way ANOVA. For example, a factor of three levels by stt;dat((2,3,4,5),(3,4,5,6),(4,3,5,7));f1 ||. Two factors A and B with 3 levels for A and 2 levels for B by stt;dat(((2,3,4,5),(8,9,8,8)),((3,4,5,6),(7,7,8,8)),((4,3,5,7),(9,10,11,9)));f2 ||. For two-way ANOVA, use "n" for interaction. Check stt;dat(((2,3,4,5),(8,9,8,8)),((3,4,5,6),(7,7,8,8)),((4,3,5,7),(9,10,11,9)));f2;n ||.

4. Nonparametric test Use "mw" (Mann-Whitney U-test) to determine if two samples come from the same population in a nonparametric way. Kruskal-Wallis H-test by "kw" extends Mann-Whitney U-test to more than two samples (analogue to one-way ANOVA). Use "w" (Wilcoxon signed-rank test) to test the null hypothesis that two related paired samples come from the same distribution (nonparametric version of paired-samples t-test). When using the Wilcoxon signed-rank test, input the differences of each pair in one row. For example, stt;dat((1,7,9,7,12,9,12,9,6,12),(3,7,7,8,9,10,8,7,6,13));mw || stt;dat((2,3,4,5),(3,4,5,6),(4,3,5,7));kw || stt;dat(-2,-1,2,-1,-2);w ||.

5. Test of independence Use "cnt" (contingency table) for test of independence. For example, check the observed counts in a
3 × 2 table by stt;dat((32,43,25),(27,30,30));cnt ||.

V. Generalized Linear Models

Use "ols" for both simple or multiple linear regression. Note that the response variable is listed in the last row and explanatory variable(s) are listed before the response variable. For example, stt;dat((2,3,6,8,1,3,5,9,2),(8,10,9,12,4,6,8,13,3));ols || stt;dat((2,3,6,8,1,3,5,9,2),(0,0,0,0,1,1,1,1,1),(8,10,9,12,4,6,8,13,3));ols ||. For predicted values and residuals, use "ols;pr"; for the predicted value given a new observation, use "ols;p(a)" or "ols;p(a,b)". For example, stt;dat((2,3,6,8,1,3,5,9,2),(8,10,9,12,4,6,8,13,3));ols;pr || stt;dat((2,3,6,8,1,3,5,9,2),(8,10,9,12,4,6,8,13,3));ols;p(4) || stt;dat((2,3,6,8,1,3,5,9,2),(0,0,0,0,1,1,1,1,1),(8,10,9,12,4,6,8,13,3));ols;pr || stt;dat((2,3,6,8,1,3,5,9,2),(0,0,0,0,1,1,1,1,1),(8,10,9,12,4,6,8,13,3));ols;p(4,1) ||.

Use "glm;b" for logistic regression, where "b" indicates the binomial family and logit link function. Note that the response variable, which is binary, must be listed in the last row in the data set. For example, stt;dat((3.78,2.44,2.09,0.14,1.72,1.65,4.92,4.37,4.96,4.52,3.69,5.88),(0,0,0,0,0,0,1,1,1,1,1,1));glm;b ||.

Use "pr;n" (polynomial regression) for polynomial regression, where "n" is a positive integer for the order of degree. For example, stt;dat((3,2,3,1,7,5,2,9),(12,5,8,2,38,28,5,90));pr;2 ||. Check the scatter plot by plt;dat((3,2,3,1,7,5,2,9),(12,5,8,2,38,28,5,90));pt ||.

Use "log" for log linear models. For three-way tables, add model generators such as "12", "23", "13", "12+23", "12+13", "13+23", "12+13+23" to the end. Note that the row and column indicators are listed as explanatory variables before the counts, which are listed in the last row as response variable. For example; stt;dat((1,1,1,1,0,0,0,0),(0,1,0,1,0,1,0,1),(0,1,0,1,0,1,0,1),(35,59,47,112,42,77,26,76));log || stt;dat((1,1,1,1,0,0,0,0),(0,1,0,1,0,1,0,1),(0,1,0,1,0,1,0,1),(35,59,47,112,42,77,26,76));log;12+23 ||.

VI. Random Process

Use "dmc" (discrete Markov chain) to create a finite discrete time-homogeneous Markov chain. The first row "(1,2,3)" in the example "stt;dmc((1,2,3),[(0.2,0.3,0.5),(0.3,0.5,0.2),(0.5,0.3,0.2)]);tp" lists the discrete state space, and the second portion is a 3 × 3 (constant) one-step transition matrix. Note that each row of the transition matrix must sum to 1. The keyword "tp" (transition probability) displays the state space and the transition matrix.

Use keyword "pn;n" for n-step transition matrix, where "n" is a positive integer for "n-step". Use "ld" for limiting distribution, "fm" for fundamental matrix, "cc" for communication class, "ap" for absorbing probabilities, and "dc" for decomposing transition matrix into submatrices with special properties. Use "md;n;(p1, ..., pk)" for n-step marginal distribution with initial distribution "(p1, ..., pk)" that sums to 1. For example, stt;dmc((0,1,2,3,4),[(1/2,1/2,0,0,0),(2/5,1/5,2/5,0,0), (0,0,1,0,0),(0,0,1/2,1/2,0),(1/2,0,0,0,1/2)]);dc || stt;dmc((0,1,2,3,4),[(1/2,1/2,0,0,0),(2/5,1/5,2/5,0,0), (0,0,1,0,0),(0,0,1/2,1/2,0),(1/2,0,0,0,1/2)]);tp || stt;dmc((1,2,3),[(0.2,0.3,0.5),(0.3,0.5,0.2),(0.5,0.3,0.2)]);pn;7 || stt;dmc((1,2,3),[(0.2,0.3,0.5),(0.3,0.5,0.2),(0.5,0.3,0.2)]);ld || stt;dmc((0,1,2,3,4),[(1/2,1/2,0,0,0),(2/5,1/5,2/5,0,0), (0,0,1,0,0),(0,0,1/2,1/2,0),(1/2,0,0,0,1/2)]);cc || stt;dmc((0,1,2,3,4),[(1/2,1/2,0,0,0),(2/5,1/5,2/5,0,0), (0,0,1,0,0),(0,0,1/2,1/2,0),(1/2,0,0,0,1/2)]);fm || stt;dmc((1,2,3),[(0.2,0.3,0.5),(0.3,0.5,0.2),(0.5,0.3,0.2)]);md;7;(1,0,0) ||.

Use "cmc" for continuous-time Markov chain. For example, stt;cmc((0,1),((-1,1),(1,-1)));P(x(1.96),0|x(0.78),1) ||.

Compute transition probabilities The expression "P(x(3),1|x(0),1)" represents the probability that the process starts 1 at time 0 and moves to 1 at time 3. Similarly, "P(x(5)>=1|x(2)>1)" means the probability that the process takes values greater than 1 at time 2 and takes values greater or equal to 1 at time 5. For example, stt;dmc((1,2,3),[(0.2,0.3,0.5),(0.3,0.5,0.2),(0.5,0.3,0.2)]);P(x(3),1|x(0),1) || stt;dmc((1,2,3),[(0.2,0.3,0.5),(0.3,0.5,0.2),(0.5,0.3,0.2)]);P(x(3),1|x(0)>=2) || stt;dmc((1,2,3),[(0.2,0.3,0.5),(0.3,0.5,0.2),(0.5,0.3,0.2)]);P(x(6)!=1|x(2)<=2) || stt;cmc((1,2,3),((-6,3,3),(4,-12,8),(15,3,-18)));P(x(3.12),3|x(1)>1) ||.