Title
mm_finvert() -- Numerical inversion of univariate function
Syntax
real vector mm_finvert(y, f, df [, x0, tol, itr])
real vector mm_finvert(y, f, lo, up [, tol, itr])
where
y: real vector containing the values for which function f be inverted
f: pointer scalar containing address of function to be inverted; usually this is coded &funcname()
df: pointer scalar containing address of function providing first derivative of function f (Newton-Raphson method only)
x0: real vector containing initial guess (Newton-Raphson method only); if x0 omitted, y is used as initial guess
lo: real vector containing lower endpoint of the search interval (Brent's method only)
up: real vector containing upper endpoint of the search interval (Brent's method only)
tol: real scalar specifying acceptable tolerance for the estimate (default is tol=0 to find the solution as accurate as possible)
itr: real scalar specifying the maximum number of iterations (default is itr=1000)
Description
mm_finvert() numerically inverts function f for the outcomes y, that is, mm_finvert() returns an approximation for x given y, where y = f(x).
Two methods are available:
mm_finvert(y, f, df [, x0, ...]): If function df, providing the first derivative of f is given, then the Newton-Raphson method implemented in mm_nrroot() is applied. x0 may be used to specify an initial guess for x. Specify x0 as a scalar to use the same initial guess for x with all outcomes in y. If x0 is omitted, y is used as the initial guess.
mm_finvert(y, f, lo, up [, ...]): If derivatives are not provided, Brent's method implemented in mm_root() is applied. With this method, lo and up, the lower and and upper endpoints of the search interval, have to be specified. Specify lo and up as scalars to use the same search interval for x with all outcomes in y.
With both methods, tol sets the acceptable tolerance for the accuracy of the approximation of x. See help for mm_nrroot() and mm_root() for details. Furthermore, itr sets the maximum number of iterations. mm_finvert() aborts with error, if convergence is not reach within itr iterations for any element in y.
mm_finvert() is loosely based on suggestions made by Jeff Pitblado on Statalist (see http://statacorp.com/statalist/archive/2005-09/msg00837.html).
Remarks
mm_finvert() may be useful, for example, generate random draws from a given distribution function. If the first derivative of the distribution function, i.e. the density function, is known, the approximation of the inversion can be computed using the Newton-Raphson method. For example, normally distributed random data could be produced as follows:
: function fx(x) return(normal(x))
: function dfx(x) return(normalden(x))
: y = uniform(5,1) : mm_finvert(y, &fx(), &dfx()) 1 +----------------+ 1 | 1.146369867 | 2 | -.3431823749 | 3 | -1.163062788 | 4 | -.9168871812 | 5 | -.3915779721 | +----------------+
If a function providing the first derivative is not available, the inversion can be computed using Brent's zero finding method. Brent's method works very well in a variety of contexts. However, it is usually slower than the Newton-Raphson method and it requires the user to specify a range within which the solution is to be sought for. Example:
: mm_finvert(y, &fx(),-5, 5) 1 +----------------+ 1 | 1.146369867 | 2 | -.3431823749 | 3 | -1.163062788 | 4 | -.9168871812 | 5 | -.3915779721 | +----------------+
The above examples have only illustrative purpose. More accurate and much more efficient would be to use the built-in function invnormal() in this context (see help for normal()), that is:
: invnormal(y) 1 +----------------+ 1 | 1.146369867 | 2 | -.3431823749 | 3 | -1.163062788 | 4 | -.9168871812 | 5 | -.3915779721 | +----------------+
Conformability
mm_finvert(y, f, df, x0, tol, itr): y: n x 1 or 1 x n f: 1 x 1 df: 1 x 1 x0: n x 1 or 1 x n or 1 x 1 tol: 1 x 1 itr: 1 x 1 result: n x 1 or 1 x n
mm_finvert(y, f, lo, up, tol, itr): y: n x 1 or 1 x n f: 1 x 1 lo: n x 1 or 1 x n or 1 x 1 up: n x 1 or 1 x n or 1 x 1 tol: 1 x 1 itr: 1 x 1 result: n x 1 or 1 x n
Diagnostics
None.
Source code
mm_finvert.mata
Author
Ben Jann, ETH Zurich, jann@soz.gess.ethz.ch
Also see