## Tag: Numerical Analysis

0
Create an interpolating polynomial in SAS

While reviewing a book on numerical analysis, I was reminded of a classic interpolation problem. Suppose you have n pairs of points in the plane: (x1,y1), (x2,y2), ..., (xn,yn), where the first coordinates are distinct. Then you can construct a unique polynomial of degree (at most) n-1 that passes through

0
Distributions with specified skewness and kurtosis

A SAS programmer wanted to simulate samples from a family of Beta(a,b) distributions for a simulation study. (Recall that a Beta random variable is bounded with values in the range [0,1].) She wanted to choose the parameters such that the skewness and kurtosis of the distributions varied over range of

Analytics
0
The geometry of Jacobi's method

A colleague remarked that my recent article about using Jacobi's iterative method for solving a linear system of equations "seems like magic." Specifically, it seems like magic that you can solve a certain class of linear systems by using only matrix multiplication. For any initial guess, the iteration converges to

0
Implement Jacobi's method in SAS

In a first course in numerical analysis, students often encounter a simple iterative method for solving a linear system of equations, known as Jacobi's method (or Jacobi's iterative method). Although Jacobi's method is not used much in practice, it is introduced because it is easy to explain, easy to implement,

0
The distribution of the sample median for normal data

A previous article shows how to use Monte Carlo simulation to approximate the sampling distribution of the sample mean and sample median. When x ~ N(0,1) are normal data, the sample mean is also normal, and there are simple formulas for the expected value and the standard error of the

0
On the evaluation of the function exp(x) - 1

Recently I wrote about numerical analysis problem: the accurate computation of log(1+x) when x is close to 0. A naive computation of log(1+x) loses accuracy if you call the LOG function, which is why the SAS language provides the built-in LOG1PX for this computation. In addition, I showed that you

0
Approximate functions by using Taylor series and rational functions

SAS supports a special function for the accurate evaluation of log(1+x) when x is near 0. The LOG1PX function is useful because a naive computation of log(1+x) loses accuracy when x is near 0. This article demonstrates two general approximation techniques that are often used in numerical analysis: the Taylor

0
Avoid domain errors by using Taylor series

The other day I was trying to numerically integrate the function f(x) = sin(x)/x on the domain [0,∞). The graph of this function is shown to the right. In SAS, you can use the QUAD subroutine in SAS IML software to perform numerical integration. Some numerical integrators have difficulty computing

0
The exponential of a matrix

In SAS, you can approximate the exponential of a matrix by using the EXPMATRIX function in SAS IML software. This article discusses the exponential of a matrix: what it is, how to compute it, why it is useful, and why you should think of it as a linear map that

0
How Newton calculated pi to 16 decimal places

Happy Pi Day! Every year on March 14th (written 3/14 in the US), people in the mathematical sciences celebrate "all things pi-related" because 3.14 is the three-decimal approximation to π ≈ 3.14159265358979.... Modern computer methods and algorithms enable us to calculate 100 trillion digits of π. However, I think it

0
The distribution of the difference between two beta random variables

A SAS programmer wanted to compute the distribution of X-Y, where X and Y are two beta-distributed random variables. Pham-Gia and Turkkan (1993) derive a formula for the PDF of this distribution. Unfortunately, the PDF involves evaluating a two-dimensional generalized hypergeometric function, which is not available in all programming languages.

0
Compute hypergeometric functions in SAS

In applied mathematics, there is a large collection of "special functions." These function appera in certain applications, often as the solution to a differential equation, but also in the definition of probability distributions. For example, I have written about Bessel functions, the complete and incomplete gamma function, and the complete

0
A test for monotonic sequences and functions

Monotonic transformations occur frequently in math and statistics. Analysts use monotonic transformations to transform variable values, with Tukey's ladder of transformations and the Box-Cox transformations being familiar examples. Monotonic distributions figure prominently in probability theory because the cumulative distribution is a monotonic increasing function. For a continuous distribution that is

0
Tips for computing right-tail probabilities and quantiles

A colleague was struggling to compute a right-tail probability for a distribution. Recall that the cumulative distribution function (CDF) is defined as a left-tail probability. For a continuous random variable, X, with density function f, the CDF at the value x is     F(x) = Pr(X ≤ x) = ∫

0
The derivative of a quantile function

Recently, I needed to solve an optimization problem in which the objective function included a term that involved the quantile function (inverse CDF) of the t distribution, which is shown to the right for DF=5 degrees of freedom. I casually remarked to my colleague that the optimizer would have to

0
Two perspectives on numerical integration

To a numerical analyst, numerical integration has two meanings. Numerical integration often means solving a definite integral such as \$\$int_{a}^b f(s) , ds\$\$. Numerical integration is also called quadrature because it computes areas. The other meaning applies to solving ordinary differential equations (ODEs). My article about new methods for solving

0
New methods for solving differential equations in SAS

Many discussions and articles about SAS Viya emphasize its ability to handle Big Data, perform parallel processing, and integrate open-source languages. These are important issues for some SAS customers. But for customers who program in SAS and do not have Big Data, SAS Viya is attractive because it is the

0
A differential equation for quantiles of a distribution

The graph to the right is the quantile function for the standard normal distribution, which is sometimes called the probit function. Given any probability, p, the quantile function gives the value, x, such that the area under the normal density curve to the left of x is exactly p. This

0
Five constants every statistical programmer should know

Statistical programmers need to access numerical constants that help us to write robust and accurate programs. Specifically, it is necessary to know when it is safe to perform numerical operations such as raising a number to a power without exceeding the largest number that is representable in finite-precision arithmetic. This

0
Least-squares optimization and the Gauss-Newton method

A previous article showed how to use SAS to compute finite-difference derivatives of smooth vector-valued multivariate functions. The article uses the NLPFDD subroutine in SAS/IML to compute the finite-difference derivatives. The article states that the third output argument of the NLPFDD subroutine "contains the matrix product J`*J, where J is

0
Finite-difference derivatives of vector-valued functions

I previously showed how to use SAS to compute finite-difference derivatives for smooth scalar-valued functions of several variables. You can use the NLPFDD subroutine in SAS/IML software to approximate the gradient vector (first derivatives) and the Hessian matrix (second derivatives). The computation uses finite-difference derivatives to approximate the derivatives. The

0
Finite-difference derivatives in SAS

Many applications in mathematics and statistics require the numerical computation of the derivatives of smooth multivariate functions. For simple algebraic and trigonometric functions, you often can write down expressions for the first and second partial derivatives. However, for complicated functions, the formulas can get unwieldy (and some applications do not

0
Test for a symmetric matrix

It is important to be able to detect whether a numerical matrix is symmetric. Some operations in linear algebra require symmetric matrices. Sometimes, you can use special algorithms to factor a symmetric matrix. In both cases, you need to test a matrix for symmetry. A symmetric matrix must be square.

0
Implement a line search algorithm in SAS

Recently, I needed to implement a line search algorithm in SAS. The line search is illustrated by the figure at the right. You start with a point, p, in d-dimensional space and a direction vector, v. (In the figure, d=2, but in general d > 1.) The goal is to

0
The ERF and ERFC functions for statisticians

Recently, a SAS programmer commented about one of my blog posts. He said that he had found an alternative answer on another website. Whereas my answer was formulated in terms of the normal cumulative distribution function (CDF), the other answer used the ERF function. This article shows the relationship between

0
Double integrals by using Monte Carlo methods

As mentioned in my article about Monte Carlo estimate of (one-dimensional) integrals, one of the advantages of Monte Carlo integration is that you can perform multivariate integrals on complicated regions. This article demonstrates how to use SAS to obtain a Monte Carlo estimate of a double integral over rectangular and

0
Sample size for the Monte Carlo estimate of an integral

A previous article shows how to use Monte Carlo simulation to estimate a one-dimensional integral on a finite interval. A larger random sample will (on average) result in an estimate that is closer to the true value of the integral than a smaller sample. This article shows how you can

0
Estimate an integral by using Monte Carlo simulation

Numerical integration is important in many areas of applied mathematics and statistics. For one-dimensional integrals on the interval (a, b), SAS software provides two important tools for numerical integration: For common univariate probability distributions, you can use the CDF function to integrate the density, thus obtaining the probability that a