R/fit_gaussian_2D.R
fit_gaussian_2D.Rd
Determine the bestfit parameters for a specific 2DGaussian 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 2DGaussian to
the supplied data. Each method uses (slightly) different sets of
parameters. Note that for a small (but nontrivial) 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 2DGaussian 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 (xaxis peak location), Y_peak (yaxis peak location),
X_sig (spread along xaxis), and Y_sig (spread along yaxis).
A more generic method (and the default) is method = "elliptical"
.
This allows the fitted 2DGaussian 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 xaxis in the
clockwise direction), X_peak (xaxis peak location), Y_peak (yaxis peak
location), a (width of Gaussian along xaxis), and b (width of Gaussian
along yaxis).
A third method is method = "elliptical_log"
. This is a further
special case in which log2transformed data may be used. See Priebe et al.
2003 for more details. Parameters from this model include: Amp (amplitude),
Q (orientation parameter), X_peak (xaxis peak location), Y_peak (yaxis
peak location), X_sig (spread along xaxis), and Y_sig (spread along
yaxis).
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 bestfit 2DGaussian (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 xaxis in the clockwise direction. In contrast, the
method = "elliptical_log"
procedure uses a Q parameter to determine
the orientation of the 2DGaussian. Setting constrain_orientation =
0
will result in a diagonallyoriented 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):565061. doi: 10.1523/JNEUROSCI.231305650.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) }