fit_gaussian_2D: Determine the best-fit parameters for a specific 2D Gaussian...

Description Usage Arguments Details Value Author(s) References Examples

View source: R/fit_gaussian_2D.R

Description

Determine the best-fit parameters for a specific 2D Gaussian model

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
fit_gaussian_2D(
  data,
  method = "elliptical",
  constrain_amplitude = FALSE,
  constrain_orientation = "unconstrained",
  user_init = NULL,
  maxiter = 1000,
  verbose = FALSE,
  print_initial_params = FALSE,
  ...
)

Arguments

data

A data.frame that contains the raw data (generally rectilinearly gridded data, but this is not a strict requirement). Columns must be named "X_values", "Y_values" and "response".

method

Choice of "elliptical", "elliptical_log", or "circular". Determine which specific implementation of 2D Gaussian to use. See Details for more.

constrain_amplitude

Default FALSE; should the amplitude of the Gaussian be set to the maximum value of the "response" variable (TRUE), or should the amplitude fitted by stats::nls() (FALSE)?

constrain_orientation

If using "elliptical" or method = "elliptical_log", should the orientation of the Gaussian be unconstrained (i.e. the best-fit orientation is returned) or should it be pre-set by the user? See Details for more. Defaults to "unconstrained".

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 stats::nls.control() for more details.

verbose

TRUE or FALSE; should the trace of the iteration be printed? See the trace argument of stats::nls() for more detail.

print_initial_params

TRUE or FALSE; should the set of initial parameters supplied to stats::nls() be printed to the console? Set to FALSE by default to avoid confusion with the fitted parameters attained after using stats::nls().

...

Additional arguments passed to stats::nls.control()

Details

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 ....

Value

A list with the components:

Author(s)

Vikram B. Baliga

References

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.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
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)
}

Example output



gaussplotR documentation built on Jan. 13, 2021, 3:55 p.m.