R/fit_gaussian_2D.R
fit_gaussian_2D.RdDetermine the best-fit parameters for a specific 2D-Gaussian model
fit_gaussian_2D( data, method = "elliptical", constrain_amplitude = FALSE, constrain_orientation = "unconstrained", user_init = NULL, maxiter = 1000, verbose = FALSE, print_initial_params = FALSE, ... )
| data | A data.frame that contains the raw data (generally rectilinearly
gridded data, but this is not a strict requirement). Columns must be named
|
|---|---|
| method | Choice of |
| constrain_amplitude | Default FALSE; should the amplitude of the
Gaussian be set to the maximum value of the |
| constrain_orientation | If using |
| user_init | Default NULL; if desired, the user can supply initial values for the parameters of the chosen model. See Details for more. |
| maxiter | Default 1000. A positive integer specifying the maximum number
of iterations allowed. See |
| verbose | TRUE or FALSE; should the trace of the iteration be printed?
See the |
| print_initial_params | TRUE or FALSE; should the set of initial
parameters supplied to |
| ... | Additional arguments passed to |
A list with the components:
"coefs" A data.frame of fitted model parameters.
"model" The model object, fitted by stats::nls().
"model_error_stats" A data.frame detailing the rss, rmse, deviance, AIC, R2 and adjusted R2 of the fitted model.
"fit_method" A character vector that indicates which method and orientation strategy was used by this function.
stats::nls() is used to fit parameters for a 2D-Gaussian to
the supplied data. Each method uses (slightly) different sets of
parameters. Note that for a small (but non-trivial) proportion of data
sets, nonlinear least squares may fail due to singularities or other
issues. Most often, this occurs because of the starting parameters that are
fed in. By default, this function attempts to set default parameters by
making an educated guess about the major aspects of the supplied data.
Should this strategy fail, the user can make use of the user_init
argument to supply an alternate set of starting values.
The simplest method is method = "circular". Here, the 2D-Gaussian is
constrained to have a roughly circular shape (i.e. spread in X- and Y- are
roughly equal). If this method is used, the fitted parameters are: Amp
(amplitude), X_peak (x-axis peak location), Y_peak (y-axis peak location),
X_sig (spread along x-axis), and Y_sig (spread along y-axis).
A more generic method (and the default) is method = "elliptical".
This allows the fitted 2D-Gaussian to take a more ellipsoid shape (but note
that method = "circular" can be considered a special case of this).
If this method is used, the fitted parameters are: A_o (a constant term),
Amp (amplitude), theta (rotation, in radians, from the x-axis in the
clockwise direction), X_peak (x-axis peak location), Y_peak (y-axis peak
location), a (width of Gaussian along x-axis), and b (width of Gaussian
along y-axis).
A third method is method = "elliptical_log". This is a further
special case in which log2-transformed data may be used. See Priebe et al.
2003 for more details. Parameters from this model include: Amp (amplitude),
Q (orientation parameter), X_peak (x-axis peak location), Y_peak (y-axis
peak location), X_sig (spread along x-axis), and Y_sig (spread along
y-axis).
If using either method = "elliptical" or method =
"elliptical_log", the "constrain_orientation" argument can be used
to specify how the orientation is set. In most cases, the user should use
the default "unconstrained" setting for this argument. Doing so will
provide the best-fit 2D-Gaussian (assuming that the solution yielded by
stats::nls() converges on the global optimum).
Setting constrain_orientation to a numeric (e.g.
constrain_orientation = pi/2) will force the orientation of the
Gaussian to the specified value. Note that this is handled differently by
method = "elliptical" vs method = "elliptical_log". In
method = "elliptical", the theta parameter dictates the rotation, in
radians, from the x-axis in the clockwise direction. In contrast, the
method = "elliptical_log" procedure uses a Q parameter to determine
the orientation of the 2D-Gaussian. Setting constrain_orientation =
0 will result in a diagonally-oriented Gaussian, whereas setting
constrain_orientation = -1 will result in horizontal orientation.
See Priebe et al. 2003 for more details.
The user_init argument can also be used to supply a vector of
initial values for the A, Q, X_peak, Y_peak, X_var, and Y_var parameters.
If the user chooses to make use of user_init, then a vector
containing all parameters must be supplied in a particular order.
Additional arguments to the control argument in stats::nls()
can be supplied via ....
Priebe NJ, Cassanello CR, Lisberger SG. The neural representation of speed in macaque area MT/V5. J Neurosci. 2003 Jul 2;23(13):5650-61. doi: 10.1523/JNEUROSCI.23-13-05650.2003.
Vikram B. Baliga
if (interactive()) { ## Load the sample data set data(gaussplot_sample_data) ## The raw data we'd like to use are in columns 1:3 samp_dat <- gaussplot_sample_data[,1:3] #### Example 1: Unconstrained elliptical #### ## This fits an unconstrained elliptical by default gauss_fit <- fit_gaussian_2D(samp_dat) ## Generate a grid of x- and y- values on which to predict grid <- expand.grid(X_values = seq(from = -5, to = 0, by = 0.1), Y_values = seq(from = -1, to = 4, by = 0.1)) ## Predict the values using predict_gaussian_2D gauss_data <- predict_gaussian_2D( fit_object = gauss_fit, X_values = grid$X_values, Y_values = grid$Y_values, ) ## Plot via ggplot2 and metR library(ggplot2); library(metR) ggplot_gaussian_2D(gauss_data) ## Produce a 3D plot via rgl rgl_gaussian_2D(gauss_data) #### Example 2: Constrained elliptical_log #### ## This fits a constrained elliptical, as in Priebe et al. 2003 gauss_fit <- fit_gaussian_2D( samp_dat, method = "elliptical_log", constrain_orientation = -1 ) ## Generate a grid of x- and y- values on which to predict grid <- expand.grid(X_values = seq(from = -5, to = 0, by = 0.1), Y_values = seq(from = -1, to = 4, by = 0.1)) ## Predict the values using predict_gaussian_2D gauss_data <- predict_gaussian_2D( fit_object = gauss_fit, X_values = grid$X_values, Y_values = grid$Y_values, ) ## Plot via ggplot2 and metR ggplot_gaussian_2D(gauss_data) ## Produce a 3D plot via rgl rgl_gaussian_2D(gauss_data) }