MLAB: An Advanced System for Mathematical and Statistical Modeling =================================================================== MLAB, (for Modeling LABoratory), is a program for interactive mathematical and statistical modeling. MLAB was originally developed at the National Institutes of Health. It includes curve-fitting, differential equations, statistics and graphics as some of its major capabilities. (Two independent reviews of MLAB may be found in IEEE Spectrum, August 1994, page 15, and Notices of the American Mathematical Society, February 1993, page 152.) MLAB provides more than thirty command types and more than four hundred built-in functions from the areas of elementary mathematics, transcendental functions, probability and statistics, linear algebra, optimization, cluster analysis, combinatorics, numeric input/output, and graphics. The usual low-level functions, e.g., sine, cosine, log, etc., are present, as well as functions performing more complex analyses, such as singular value decomposition, discrete Fourier transforms, solution of differential equation systems, and constrained non-linear optimization, among many others. A substantial collection of statistically-oriented functions, such as most common distribution functions and their inverses, are included. MLAB is an ideal tool for solving simulation and modeling problems such as chemical kinetics, pharmocological compartmental models, multiple site ligand binding, neurophysiological modeling, and ultracentrifuge models, to name just a few. Using the built-in functions in MLAB, the user may construct user-defined functions, which may be algebraic, differential equation-defined, implicit, recursive, and/or defined piecewise. Functions may be used in many ways: direct evaluation, symbolic derivatives may be constructed, and roots of nonlinear equations may be found. - MLAB can solve systems of differential equations. These differential equations may be of any order, stiff or non-stiff, and may contain delays. MLAB uses an algorithm which automatically selects both method and stepsize to minimize both truncation error and computer time. However, control variables are provided to permit the user to specify the method of numerical integration, the error tolerance, and the method for dealing with singularities. - A major capability of MLAB is constrained, non-linear optimization. Data may be input which is to be fit by one or more model functions which contain one or more parameters. Any type of function, including differential equation-defined functions, may be used as a model. Once the parameters are initialized, The fit command will adjust the parameters so as to minimize the error between the data and the model. The parameters may appear nonlinearly in the model, and be subject to multiple linear equality and inequality constraints. Several control variables are provided to permit the user to fine-tune the algorithm, although the defaults are quite robust. - MLAB provides powerful and flexible 2D and 3D graphics facilities. With a single simple command, a graph may be created with axes. A few more commands serve to insert titles and axis labels. The user has complete control of color, line-types and point-types, fonts and the position of every aspect of the plot. Many types of plots are available, including linear, logarithmic, and polar coordinates. Contour plots for functions of two variables are provided. Screen plots may be saved as files suitable for printing on HP or Postscript laser printers, at the full resolution of the output device. - MLAB is also a programming language. It provides facilities for branching, looping, input/output, and interruption of execution. Users may write MLAB program files containing multiple commands, including interactive commands, and calls to other such files. Analyses may thus be programmed much more compactly than with conventional programming languages because of the very high level of the individual commands. - MLAB is extraordinarily polished. Careful attention has been paid to significant details so that MLAB operates flawlessly, just as you think it should behave from a mathematical viewpoint. For example, the MLAB matrix sorting operator, sort(m,j), returns a copy of the matrix m with its rows sorted in ascending order. The sort is stable, so that rows with duplicate values in column j appear in the output in the same relative order as they appeared in the input. If j is negative, the sorting is done so that column |j| appears in descending order. If j is not given, column 1 is assumed, so sort(m) has a natural meaning. Finally if j is 0, the output will be a copy of m with its rows stably sorted in lexicographic order. Another very important example of attention to detail in MLAB is that MLAB correctly handles floating point overflows and underflows, supplying whichever of MAXPOS, the largest computational value, or MAXNEG, the algebraically least computational value, is indicated as the result of an overflow (with an optional warning). It supplies 0 as the result of a hard underflow. Often-appropriate values are supplied for zero-divides as well. Since overflows may well arise during curve-fitting, this corrective behavior is crucial to allowing the curve-fit to continue successfully. MLAB runs on IBM PC's and compatibles using the DOS operating system. It requires a math co-processor, VGA or EGA color graphics, and a laser printer for hard copy. UNIX versions for Sun Sparc workstations, NeXT machines, and SGI workstations are available. A Macintosh version is also available. A Windows version is under development. MLAB may be purchased from: Civilized Software, Inc. 7735 Old Georgetown Rd. #410 Bethesda, MD 20814 Tel: (301) 652-4714 Fax: (301) 656-1069 Email: csi%sava@cs.umd.edu A three month trial version is available for $100. The full system for DOS with technical support and updates for one year is priced at $1495, the Intel-NeXT system is $1495, the Motorala NeXT system is $1995, the Machintosh version is $1995, UNIX versions ( SUN, SGI) are $2995, and a streamlined DOS-version is available for $995. There is a discount of min(10*(n€1),60)% for a purchase of n copies. And 1-year lifetime copies are available for $100 apiece to educational purchasers who have at least one full MLAB copy. __________________________________________________________________ An MLAB Example: The Dynamics of Radioactive Decay Radioactive decay is a process which occurs at a rate characteristic of the decaying isotope and proportional to the amount of the isotope. The decay process converts the original isotope into a new species called a daughter, which, itself, may also be radioactive and undergo further decay, again at a characteristic rate proportional to the amount of the second isotope, which we may call the granddaughter species. The granddaughter may also be radioactive, but have a non-radioactive decay product. From observation of the time-dependence of the radioactivity of the grand-daughter, we wish to estimate the decay rates of all the three species, and hence, to identify the three isotopes which are involved. This type of identification problem was important, for example, in determining the atomic composition of the star involved in the recently observed supernova. Denote the amounts of the three isotopes by X, Y, and Z, respectively. Denote the initial amounts of the isotopes by X0, Y0, and Z0; and denote the corresponding rates of decay by K, L, and M, which are positive numbers. The differential equations of the radioactive decay process described above are: dX/dt = €K*X dY/dt = K*X€L*Y dZ/dt = L*Y€M*Z We suppose that only the decay of isotope Z produces radiation which can be detected, and that measurements of radiation from decays of isotope Z (which is proportional to the amount of Z) are made at a series of times Ti, i = 1,2,...,200. These measurements are read into a 200 row, two column matrix called DATA, column 1 of which contains the times, and column 2 of which contains the measured values. We suppose that the initial conditions are known, e.g. X0 = 5000, Y0 = 0, and Z0 = 0. Initial guesses are required for K, L, and M, from which the process of estimating K,L, and M will proceed. We take these guesses to be K=1, L=1, and M=1. In MLAB, this problem may be formulated as follows: * FCT X'T (T) = €K*X * FCT Y'T (T) = K*X€L*Y * FCT Z'T (T) = L*Y€M*Z * INIT X(0) = 5000 * INIT Y(0) = 0 * INIT Z(0) = 0 * K = 1; L = .1; M = .01; * CONSTRAINTS Q = {K > 0, L > 0, M > 0} * DATA = READ(DECAY,200,2) Now we may estimate K, L, and M using the following MLAB FIT statement. * FIT (K,L,M),Z TO DATA, CONSTRAINTS Q The output of the fit statement is the table shown below, containing the least squares estimates, along with standard deviations which help the user to determine the quality of the fit. final parameter values value error dependency parameter 0.1133576575 0.9382481192 0.9996464477 K 0.1171530638 0.8648505987 0.9996463229 L 0.01384391436 0.0005651783173 0.07687919473 M 17 iterations CONVERGED best weighted sum of squares = 1.315000e+06 weighted root mean square error = 2.564176e+02 weighted deviation fraction = 9.784715e-02 R squared = 9.595270e-01 no active constraints A graph of the data and the fitted curve is a useful measure of the quality of the fit. This is easily done in MLAB. * DRAW DATA LINETYPE NONE, PT CIRCLE, PTSIZE 0.01, COLOR RED * DRAW POINTS(Z,0:200!160) COLOR BLUE * IMAGE COLOR GREY * BOTTOM TITLE "time" * LEFT TITLE "Z-isotope amount" * TOP TITLE "Cascaded exponential decay curve fit" * VIEW MLAB is unique in its power and ease of use in handling differential equations in curve-fitting. Delays, implicitly-defined functions, and recursively-defined functions may all be used. The MLAB differential equation solver uses Adams' fourth order predictor-corrector method for speed and Gear's implicit predictor-corrector method for stiff equations. The derivatives needed for the elements of the Jacobian in Gear's method are computed symbolically and evaluated from the resulting formulas. Integration step size is chosen automatically and varied so as to maximize speed while maintaining specified accuracy. __________________________________________________________________