help mata mm_finvert()-------------------------------------------------------------------------------

Title

mm_finvert() -- Numerical inversion of univariate function

Syntax

real vectormm_finvert(y,f,df[,x0,tol,itr])

real vectormm_finvert(y,f,lo,up[,tol,itr])where

y:real vectorcontaining the values for which functionfbe inverted

f:pointer scalarcontaining address of function to be inverted; usually this is coded&funcname()

df:pointer scalarcontaining address of function providing first derivative of functionf(Newton-Raphson method only)

x0:real vectorcontaining initial guess (Newton-Raphson method only); ifx0omitted,yis used as initial guess

lo:real vectorcontaining lower endpoint of the search interval (Brent's method only)

up:real vectorcontaining upper endpoint of the search interval (Brent's method only)

tol:real scalarspecifying acceptable tolerance for the estimate (default istol=0 to find the solution as accurate as possible)

itr:real scalarspecifying the maximum number of iterations (default isitr=1000)

Description

mm_finvert()numerically inverts functionffor the outcomesy, that is,mm_finvert()returns an approximation forxgiveny, wherey=f(x).Two methods are available:

mm_finvert(y,f,df[,x0,...]): If functiondf, providing the first derivative offis given, then the Newton-Raphson method implemented inmm_nrroot()is applied.x0may be used to specify an initial guess forx. Specifyx0as a scalar to use the same initial guess forxwith all outcomes iny. Ifx0is omitted,yis used as the initial guess.

mm_finvert(y,f,lo,up[,...]): If derivatives are not provided, Brent's method implemented inmm_root()is applied. With this method,loandup, the lower and and upper endpoints of the search interval, have to be specified. Specifyloandupas scalars to use the same search interval forxwith all outcomes iny.With both methods,

tolsets the acceptable tolerance for the accuracy of the approximation ofx. See help formm_nrroot()andmm_root()for details. Furthermore,itrsets the maximum number of iterations.mm_finvert()aborts with error, if convergence is not reach withinitriterations for any element iny.

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 fornormal()), 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 x1 or 1x nf: 1x1df: 1x1x0:n x1 or 1x nor 1x1tol: 1x1itr: 1x1result:n x1 or 1x n

mm_finvert(y,f,lo,up,tol,itr):y:n x1 or 1x nf: 1x1lo:n x1 or 1x nor 1x1up:n x1 or 1x nor 1x1tol: 1x1itr: 1x1result:n x1 or 1x n

DiagnosticsNone.

Source codemm_finvert.mata

AuthorBen Jann, ETH Zurich, jann@soz.gess.ethz.ch

Also see