USING MFIT This documentation is fairly brief at the moment. We hope to include further information, including real- and Fourier-space mathematical descriptions of the model components, in a future release. Running mfit/clfit/fitgui ========================= There are three interfaces to the model-fitting code. 'mfit' should be run in a terminal window, and prompts for input as it goes along. 'clfit' (for command-line fit) gets all required parameters from the command line and does not prompt for further input. Both can make plots to any specified PGPLOT device. 'fitgui' is a graphical front-end to clfit. The command-line arguments accepted by clfit can be listed by typing 'clfit' without any arguments. The prompts in mfit and the labelling of GUI elements in fitgui should be self-explanatory. If not, a little trial and error may help! Plotting ======== mfit plots the uv coverage immediately after reading the data. It also plots squared visibility, triple product amplitude, and triple product phase against the (longest) projected basline, first for the initial model, then for the best-fit model. Flagged points are shown in red. If all points of a particular data type are flagged, or there is no data of that type, the plot is skipped. After each plot, the program asks 'enter x-axis range for replot ([return] to skip)'. Enter two numbers followed by RETURN to plot the same data over a narrower x-axis range, or just RETURN to continue. clfit will generate at most one plot for each invocation. Use the '-p' argument to specify the plot type, or the '-z' argument to specify the plot type and x-axis range. Valid plot types are: uv uv coverage for squared-visibility data vis2 squared visibility against projected baseline t3amp triple product amplitude against longest projected baseline in triangle t3phi triple product phase against longest projected baseline in triangle fitgui will similarly make up to one plot each time 'Go' is clicked. The plot types are as listed above. Data files ========== mfit and clfit (and hence fitgui) will read the following types of data file: .fits OI Exchange format - see http://www.mrao.cam.ac.uk/~jsy1001/exchange/ Only the squared visibility and triple product data types are used. The file must have extension .fits [The file formats below have been used for COAST data, but won't be of much interest to the community. Details are in the COAST Data Reduction Manual] .vis Old analyse format (no errors on data points), visibility amplitudes only .nvis New analyse format (errors on data points), visibility amplitudes only .mapdat Mapping data format. Can contain (squared) visibility amplitudes and/or triple products Both .vis and .nvis files contain *signed* visibility amplitudes. The models fit by analyse are all centro-symmetric, hence have purely real visibilities, so you are supposed to assign appropriate signs (ideally using contemporaneous closure phase data) to sqrt(V^2) values when constructing a vis or nvis file. mfit ignores these signs and works with squared visibilities. You should put any closure phase measurements in explicitly (you will need to create a mapdat or OI-FITS file for this). Model file format ================= The program reads through the model file until it comes across the keyword 'source' at the begining of a line. It then reads the next 2 words as the source name. There must be 2 words present like 'alp' 'Ori' etc. The program then reads in the various components of the model. It runs through the model file until it comes across the keyword "component" at the begining of the line. It then reads the data for that particular component: N.B. to disregard a component simply replace "component" with "!component" etc - it will then be ignored as the keyword marking a component info block. Component data must follow this prescription: keyword is present on the start of the line followed by the values associated with it. Maximum number of components is set to 10. Maximum order of taylor and gauss-hermite LD types is also set to 10. Keyword Parameter information ------- --------------------- name A one word name for the component itself, e.g. 'main' or 'hot_spot' or something. shape_type One word specifying the shape of the component. Can use either 'ellipse' 'disc' or 'point' obviously corresponding to elliptical, circular or unresolved point objects. ld_type One word specifying limb darkening type of the component, 'uniform', 'taylor', 'hestroffer', 'square-root', 'gauss-hermite' or 'gaussian'. Hankel transform functions all as defined in the maths work. if 'point' is specified as the shape_type then nothing should be put here. ld_order A single integer (1-10) specifying the order of the expansion in the 'taylor' and 'gauss-hermite' LD cases. No value should be put here in other cases. position Two values follow describing the position of the component with respect to the chosen phase reference point. First number is the polar radius r of the component in mas. Second number is the polar angle theta in degrees (measured N>E) position_prior Two values giving the prior widths for the position data above (prior assumed gaussian of standard deviation = prior width) flux One value giving the flux B of the source. Visibilties are always normalised so that unit visibility is obtained on zero baselines therefore flux is a relative quantitity and has no meaning for a single component model. flux_prior One value giving prior width on the flux shape_param Up to three values giving shape info. if 'point' is the shape type nothing should be put here. if 'disc' is the shape type only one number should be put, the diameter, a, of the disc in mas. ("diameter" is defined slightly differently depending on the limb darkening model, e.g. it is FWHM for gaussian - see the maths stuff). if 'ellipse' is the shape type all three numbers are needed giving the major axis, a, in mas, the orientation, phi, of the major axis (measured N>E) in degrees and the eccentricity factor epsilon. shape_param_prior A prior for each of the numbers given above ld_param Up to twenty limb darkening parameters may be put here. For 'uniform' and 'gaussian' LD types nothing should be here. For 'hestroffer' only one paramter, the hestroffer parameter "alpha" should be put. For 'square-root' two parameters, the "alpha" and "beta" coefficients of (1-mu) and (1-root(mu)) should be put. 'gauss-hermite' and 'taylor' cases require the same number of parameters as specified in the ld_order field above. ld_param_prior Prior widths associated with each number supplied above Illegal minimisations ===================== Fitting is performed by attempting to minimise the negative logarithm of posterior probability. Parameters with non-zero priors are considered to be free parameters in the fitting process. The fitting routine will trigger an error on the following "illegal" types of minimisation: 1. Attempting to fit a one component model with position angle theta and/or position radius r as free parameter(s) 2. Attempting to fit any components that have non-free position radius r that is set to zero but have position angle theta as a free parameter (if r is zero but free to be varied then this is okay) 3. Attempting to fit a one component model with flux B as a free parameter 4. Attempting to fit any elliptical components that have non-free eccentricity set to unity but have orientation angle phi free to vary 5. Any model with no free parameters(!) 6. Attempting to fit non-centrosymmetric models to vis/nvis file types Minimisation method =================== The program attempts to minimise the unnormalised negative log posterior probability (equal to the sum of the negative log prior and negative log likelihood) of the model given the data. The data consist of any number of squared-visibilities, triple product amplitudes, and triple product phases, with their (assumed Gaussian) one-sigma errors. The minimiser used is the routine PDA_UNCMND from the Starlink PDA library (http://star-www.rl.ac.uk/star/docs/sun194.htx/sun194.html). This works better than the simplex routine from Numerical Recipes that was used in the first version of mfit. From the book, "Numerical Methods and Software" by D. Kahaner, C. Moler, S. Nash Prentice Hall, 1988: This routine uses a quasi-Newton algorithm with line search to minimize the function represented by the subroutine FCN [the negative log posterior]. At each iteration, the nonlinear function is approximated by a quadratic function derived from a Taylor series. The quadratic function is minimized to obtain a search direction, and an approximate minimum of the nonlinear function along the search direction is found using a line search. The algorithm computes an approximation to the second derivative matrix of the nonlinear function using quasi-Newton techniques. More information is available in the reference: R.B. Schnabel, J.E. Koontz, and BE.E. Weiss, A modular system of algorithms for unconstrained minimization, Report CU-CS-240-82, Comp. Sci. Dept., Univ. of Colorado at Boulder, 1982. Posterior, Goodness-of-fit ========================== At the solution position the (unnormalised) negative log posterior is given (i.e. the sum of the negative log prior and negative log likelihood without the subtraction of the negative log evidence). The chi squared value is also reported: chi sqrd = sum { ( (data-model)/(error in data) )^2 } The number of degrees of freedom (vis + triple data points - no of free parameters) is reported along with the chi sqrd value divided by this. Note that closure phases that are linear combinations of other measured closure phases are incorrectly counted as independent data points. Errors ====== The hessian matrix Hes[ij] = (d^2/(dxi*dxj))(-log posterior) is computed at the solution point. Hessian elements are estimated using a numerical differentiation method: 1. for on-diagonal elements H[ii] = (d^2/dxi^2)(-log posterior) the element is estimated to be: H[ii] = { P(x+dx) + P(x-dx) - 2P(x) }/dx^2 where x = x[i] dx = delta.x[i] 2. for off-diagonal elements H[ij] = (d^2/(dxi.dxj))(-log posterior) the element is estimated to be: H[ij] = {P(x+dx,y+dy)+P(x-dx,y-dy)-P(x+dx,y-dy)-P(x-dx,y+dy)}/{4.dx.dy} where x=x[i] y=x[j] The covariance matrix Cov[ij] is the inverse of the Hes matrix. Some problems appear to arise where diagonal elements of the matrix, corresponding to the variances in the parameters, are negative. This could possibly point to poorly chosen priors or a poor fit/irregularities in the negative log posterior surface. Where such negative values occur the program warns the user. The correlation matrix Cor[ij] = Cov[ij] / sqrt { Var[i].Var[j] } is also computed. n.b. Var[i] is defined as Cov[ii]. Problems arise in the case of negative elements on the Cov diagonals i.e. negative Var[i]'s - in this case the correlation matrix will be incomplete with mathematically illegal elements set to zero. The "hessian-based error" of a parameter is the square root of the variance in the quantity. Model selection =============== An approximation to the negative log evidence is also reported: -log(posterior) - N/2 log(2pi) + 1/2 log(det(H)) where N is the number of model parameters varied. This expression is only valid if the posterior pdf is unimodal, and the peak in the posterior can be approximated by a Gaussian. Known deficiencies ================== Other plot choices (e.g. versus hour angle - need to read OI_TARGET for this) should be added. Ignores complex visibilities from OI-FITS file. Assumes single target in OI-FITS file. $Id: documentation,v 1.8 2003/07/18 17:50:23 jsy1001 Exp $ Local variables: fill-column: 79 End: