MOdel

See definition subtopic for help on creating models.

MOdel ? List all built-in model components.

MOdel @filename Cause the model definition and parameters to be read from the file with name "filename". Although chi^2 is written to the model file it is important to refit before issueing an UNcertain command. This is to avoid potential problems with roundoff errors.

definition

A model consists of several components that are added together. For example, "MOdel CONS LINR QUAD" will add a constant term, a linear term, and a quadratic term.

For each parameter required by the "MOdel" command, you will be prompted for four numbers --- "VAL", "SIG", "PLO", and "PHI" --- as described below. For each parameter, you should enter an initial value for "VAL"; but you can usually default on the other three numbers.

"VAL": This is the actual value of the parameter. Although CURFIT will often find the the best set of parameters to model the data, it never hurts to start it with parameters near the expected best fit.

"SIG": Any value of "SIG">=0 will not affect the outcome of "Fit". After you fit the model, "SIG" will contain the one-sigma curvature errors. This number is used by the "Uncertainty" command to start a formal error determination. If the "Uncertainty" command fails to converge because the original error estimate is wrong, sometimes you can improve the convergence by adjusting "SIG" to be a better estimate before using "Uncertainty". If you set "SIG=-1", then the parameter is frozen such that CURFIT is not allowed to change the parameter value while fitting. If you set "SIG"=-IPAR, the next number ("PLO") will default to 1, such that the current parameter value is forced to equal the value of parameter IPAR. (Note: IPAR can not equal 1 or the current parameter number.) If you place a number (N) after "SIG", this will force the current parameter to be N times the specified parameter. (N defaults to 1.0.)

PLT> MOdel GAUS GAUS 1, GC: VAL( 1.00), SIG( 0.00 ), PLO( 0.00), PHI( 0.00)? ,-4,2 2, GW: VAL( 1.00), SIG( 0.00 ), PLO( 0.00), PHI( 0.00)?

3, GN: VAL( 1.00), SIG( 0.00 ), PLO( 0.00), PHI( 0.00)?

4, GC: VAL( 1.00), SIG( 0.00 ), PLO( 0.00), PHI( 0.00)? 3. 5, GW: VAL( 1.00), SIG( 0.00 ), PLO( 0.00), PHI( 0.00)?

6, GN: VAL( 1.00), SIG( 0.00 ), PLO( 0.00), PHI( 0.00)?

defines a model consisting of two gaussians, with the x values of the centers differing by a factor of 2. Although you did not enter a value for parameter 1, it will be set to the value of 6 (2 times value of parameter 4). This relation will be maintained throughout a fit.

PLO, PHI: If SIG>=0 and if PLO<PHI, the parameter value is constrained to lie in the range PLO to PHI. Note: If PLO and PHI are both the same (say equal to zero), then the parameter will not be constrained in any way. In general, if you have difficulty fitting some data, the best thing to do is to freeze some parameters near to their expected values and then fit the reduced parameter set. When a good fit has been found with the reduced set, thaw some of the parameters and refit. If this method does not work, then you may be forced to use PLO and PHI to limit certain parameters to a meaningful range.

PLT> MOdel ? will list all builtin model components.

2D_models

PLT can fit functions to two dimensional data. In general this is best done with a COD file using the Y variable. A few common, simple functions are provided below.

CGau

Circular Gaussian:

FNY=FNY+Gn*EXP(-[xs**2-ys**2]/2.)

where xs=(X-Xc)/Gw, ys=(Y-Yc)/Gw, and with integral SQRT(2*PI)*Gn*Gw*Gw.

NCga

Normalized Circular Gaussian. Similar to CGau except the norm is now the integral:

FNY=FNY+Gn*EXP(-[xs**2-ys**2]/2.)/(2.*PI*Gw*Gw)

where xs=(X-Xc)/Gw, ys=(Y-Yc)/Gw, and with integral SQRT(2*PI)*Gn*Sx*Xy.

EGau

Ellipical Gaussian:

FNY=FNY+Gn*EXP(-[xs**2-ys**2]/2.)

where xs=(X-Xc)/Sx, ys=(Y-Yc)/Sy, and with integral SQRT(2*PI)*Gn*Sx*Sy

NEga

Normalized Ellipical Gaussian. Similar to EGau except the norm is now the integral:

FNY=FNY+Gn*EXP(-[xs**2-ys**2]/2.)

where xs=(X-Xc)/Sx, ys=(Y-Yc)/Sy, and with integral Gn.

LY

Linear in Y

FNY=FNY+Ly*Y

Note the general model for a plane is

Li*(X-Xc)+Ly*(Y-Yc)+b

Expand and collect terms to get

Li*X+Ly*Y+b-Li*Xc-Ly*Yc

Which would be created with a MOdel of 'LI LY CO' and the CO would be set to an initial value of b-Li*Xc-Ly*Yc.

Note, the maximum slope of the plane would be in the direction ATAN2(Ly,Li).

CONS

Select a model with a constant component:

FNY=FNY+CO.

LINR

Select a model with a linear component:

FNY=FNY+Li*X.

QUAD

Select a model with a quadratic component:

FNY=FNY+QU*X**2.

CUBI

Select a model with a cubic component:

FNY=FNY+CU*X**3.

X4

Select a model with an x^4 component:

FNY=FNY+X4*X**4.

X5

Select a model with an x^5 component:

FNY=FNY+X5*X**5.

POWR

Select a model with a power-law component:

FNY=FNY+PN*X**IN.

SIN

Select a model with a sinusoidal component:

FNY=FNY+SN*SIN(2*PI*(X-PH)/PE).

GAUS

Select a model with a gaussian component:

FNY=FNY+GN*EXP(-Z*Z/2.),

where Z=(X-GC)/GW and with integral SQRT(2*PI)*GN*GW.

EXP

Select a model with an exponential component:

FNY=FNY+EN*EXP(-(X-EC)/EW).

AEXP

Select a model with a symmetric exponential component (exp(-|x|) for all x):

FNY=FNY+EN*EXP(-ABS(X-EC)/EW).

BURS

Select a model with a burst component (linear rise followed by an exponential decay):

FNY=FNY+0 for X<ST; FNY=FNY+BN*(X-ST)/(PT-ST) for ST<X<PT; and FNY=FNY+BN*EXP(-(X-PT)/DT) for PT<X.

SBUR

Select a model with a smooth-burst component:

FNY=FNY+BN*(T^RR)*EXP(-(X-TS)/DT),

where T=EXP(1)*(X-TS)/(RR*DT), such that SBUR = BN at the peak.

PEAR

Select a model with a Pearson-function component:

FNY=FNY+K*(F1^M1)*(F2^M2),

where F1=[1.+(X-X0)/A1] and F2=[1.-(X-X0)*M1/(A1*M2)].

WIND

Select a model with a window-function component:

FNY=FNY+LE for T1<X<T2; and FNY=FNY+0 otherwise.

KING

Select a model with a King-profile component:

FNY=FNY+S0*(1.+(X/RC)^2)^(-IN).

LORE

Select a model with a Lorentz-profile component:

FNY=FNY+LN/(1.+[ 2.*(X-LC)/LW ]^2),

where LW is the full width at half maximum (FWHM) and the integral is PI*LN*LW/2.

SPLN

Select a #-knot spline component. The number of knots defaults to 2, which generates a straight line. The spline is computed by solving equation 3.3.7 in the 1988 edition of Numerical Recipes, by Press, Flannery, Teukolsky, and Vetterling.

For unconstrained y values, the natural spline condition is imposed, which sets y''=0 at the boundaries. You may not extrapolate this function outside the interval fitted.

It is possible to impose a periodic boundary condition on the spline curve. To do this, constrain the y position of the last knot to be the same as the first. When this constraint is detected, the program automatically forces the first derivatives to match at the two boundaries. For this case, you are allowed to access the function outside the interval fitted. However, the function is assumed to be periodic, with the period given by the difference in x between the first and last knots.

For example, "MOdel SPLN 5" will generate a 5-knot spline (10 parameters). The spline can be added to other models; thus "MOdel SPLN 5 GAUS" would add a 5-knot spline to a gaussian. Hence, the spline would model the `background' and the gaussian, a `line'.

problems

It is possible for the x position of two knots to lie between two adjacent data points. This results in a local chi^2 minimum as the lower knot adjusts to fit data below it, the upper knot adjusts to fit data above it. A strong wiggle occurs between the two knots but since there are no data points there, chi^2 is not affected. In this section, two knots very close to each other will be called a collision. If collision occurs during a fit, then convergence will be very slow.

One method to greatly reduce the number of collisions is to first fit the y locations before attempting to fit the x locations. By default, the knots are evenly spaced in the x direction and are not allowed to vary. For the first "Fit" you should leave the x positions frozen, although you can move the knots (using "Newpar") to concentrate them where the function is changing rapidly. Once a reasonable set of y positions is determined you can then thaw the x positions and re-fit. You should never thaw the end points: They determine the range over which the spline is to be evaluated.

With the above recipe, collisions can still occur. The straight-forward method to separate the knots is: Use "Newpar" to re-position the two knots, freeze the x locations, and then re-fit. After this the knots will sometimes stay separated when you thaw their positions and re-fit. The trick is to force the knots far enough apart so that they will not be attracted to the local minimum, but not so far apart as to grossly distort the fit.

Sometimes two knots collide when you are trying to fit the data with too few knots. This case can be easily tested for by increasing the number of knots and re-fitting.

AKIM

Select a #-knot Akima component. An Akima component is very similar to "SPLN" in that both use a cubic function to interpolate between the knots. Akima's method does not introduce false extrema and inflection points as does the cubic spline and therefore, is far superior for data that show abrupt transitions. Details of Akima's method can be found in ``A new method of Interpolations and smooth curve fitting based on local procedures'' by Hiroshi Akima in J. of the Ass. for Computing Machinery, (1970) 17, 589-602. (An implementation is described in PPC Journal, (1985) 12, no. 10, 11-14).

Like "SPLN" two different boundary conditions are allowed. If the last y value is unconstrained, then the code uses `virtual' knots outside the boundaries to determine the function at the boundaries. The locations of the virtual knots mirror the location of the knots just inside the boundaries. If the y position of the last knot is constrained to match the y position of the first knot then a periodic boundary condition is imposed.

problems

DEMO

Call the Fortran user-defined component. The QDP/PLT User's Guide describes how to create how to write a Fortran function that can be linked in to PLT to replace the "DEMO" component.

$codfile

Call the user-defined COD (COmponent Definition) function found in "codfile.COD". Briefly, a COD function is a program written in a Forth-like computer language. To understand COD, read the documentation or the on-line help for COD. A COD file can be added to any combination of built-in components. For example, the model specified by "MOdel CONS LINR $TEST" would calculate the sum of a constant term, a linear term, and the value of the COD function contained in the file "TEST.COD".

At the present time only one COD function can be defined in a model, although this function can be referenced more than once. If you wish to combine two COD functions, you will need to write a third function that combines the first two.

COD should be used for all simple components that cannot be expressed by adding together the built-in components. Since a COD function is interpreted, it will run slower than the user-defined component. However, since COD is highly efficient and supports many mathematical functions, it is expected that the interpreter will be good enough for most purposes. For large numbers of points (> 10^4) or models that involve reading a disk file, the user is advised to write a Fortran function using the user component.

example

: NGAUS ! The file must contain a : followed by a dummy name
X       ! Push current value of X onto the stack
X       ! Push current value of X onto the stack
*       ! Multiply the top two numbers on the stack to get X*X
P1      ! Push the value of parameter 1 onto the stack
*       ! Multiply to get P1*X*X
NEG     ! Negate the number on the top of the stack (-P1*X*X)
EXP     ! Calculate EXP of -P1*X*X
P2      ! Push the value of parameter 2 onto the stack
*       ! Multiply to get P2*EXP(-P1*X*X)
;       ! The function must end with a ; character

This simple COD function ("NGAUS.COD") contains two parameters and calculates the value of "P2*EXP(-P1*X*X)". It could be written much more concisely as

: NGAUS X X * P1 * NEG EXP P2 * ;

Return to plt main page.