# Is least squares the standard method to fit a 3 parameters Gaussian function to some x and y data?

A participant in one experiment needs to decide whether a flash and a sound are simultaneous or not for many possible asynchronies between the flash and the sound (x in seconds). For each asynchrony, the flash and sound are presented 100 times. In the below graph the proportion of 'simultaneous' responses (y) is plotted as a function of the asynchrony I want to fit a 3-parameter Gaussian distribution to these data. I used least squares (below). My question is whether least squares is the most standard/orthodox form to proceed for this kind of fitting. It makes sense to use maximum likelihood estimation here? How can be done?

  x<-seq(-.3,.3,.1)
yes<-c(10,45,90,100,60,10,5) #simultaneous responses
y<-yes/100

leastSquares<-function(p) sum((y-p*dnorm(x,p,p))^2)
p<-optim(c(.001,0.1,.210),leastSquares)$par xseq<-seq(-.3,.3,.01) yseq<-p*dnorm(xseq,p,p) plot(x,y) lines(xseq,yseq)  • Fitting a Gaussian to data like these is described at stats.stackexchange.com/questions/70153/… (including R code). Least squares fits are maximum likelihood estimates for Normally-distributed residuals. As Nick Cox points out, your residuals do not have such distributions. – whuber Sep 24, 2013 at 13:31 ## 2 Answers Use logistic regression with a Gaussian link. The count of simultaneous responses for a given value of$x$, written$y(x)$, is the outcome of$n=100$independent Bernoulli trials whose chance of success is given by the Gaussian function. Letting$\theta$stand for the three parameters (unknown, to be estimated), let's write the value of that Gaussian at$x$as$\mu(x, \theta)$. Then the probability of observing$y(x)$successes, with$0 \le y(x) \le n$, is $${\Pr}_\theta(y(x)) = \binom{n}{y(x)} \mu(x,\theta)^{y(x)} \left(1-\mu(x,\theta)\right)^{n-y(x)}.$$ The likelihood of a set of independent observations arising from the same underlying Gaussian is the product of these expressions. The part of the product that varies with$\theta$is obtained from the last two factors. They will tend to be extremely small (because we are dropping the binomial coefficients), so about the only reasonable way to handle them on a computer is through their logarithms. Thus the part of the log likelihood that varies with$\theta$is $$\Lambda(\theta) = \sum_{x}\left(y(x)\log(\mu(x,\theta)) + (n-y(x))\log(1-\mu(x,\theta))\right).$$ Logistic regression with a Gaussian link maximizes this log likelihood. To make sure the Gaussian peak does not exceed$1$, we might choose to parameterize its amplitude using a function that ranges from$0$to$1$. Here is one convenient parameterization: $$\mu(x, \theta) = \mu(x, (m,s,a)) = \frac{\exp(-\frac{1}{2} \left(\frac{x-m}{s}\right)^2)}{1 + \exp(-a)}.$$ The parameter$m$is the mode of the Gaussian,$s$(a positive number) is its spread, and$a$(some real number) determines the amplitude, increasing with increasing$a$. Use a multivariate nonlinear optimization procedure suitable for smooth functions. (Avoid specifying any constraints at all by using$s^2$instead of$s$as a parameter, if necessary.) Given some vaguely reasonable estimates of the parameters, it should have no trouble finding the global optimum. As an example, here is a detailed implementation of the fitting procedure in R using data from the question. It is modified from code for a four-parameter least-squares fit of a Gaussian shown in an answer at Linear regression best polynomial (or better approach to use)?. The fit is good: the standardized residuals do not become extreme and given the small amount of data, they are reasonably close to zero. The estimated values are$\hat{m} = -0.033$,$\hat{s} = 0.127$, and$\hat{a} = 16.8$(whence the estimated amplitude is$1/(1+\exp(-\hat{a})) = 1 - 0.00000005$). y <- c(10,45,90,100,60,10,5) # Counts of successes at each x. N <- length(y) n <- rep(100, N) # Numbers of trials at each x. x <- seq(-3,3,1)/10 # # Define a Gaussian function (of three parameters {m,s,a}). # gaussian <- function(x, theta) { m <- theta; s <- theta; a <- theta; exp(-0.5*((x-m)/s)^2) / (1 + exp(-a)) } # # Compute an unnormalized log likelihood (negated). # y are the observed counts, # n are the trials, # x are the independent values, # theta is the parameter. # likelihood.log <- function(theta, y, n, x) { p <- gaussian(x, theta) -sum(y*log(p) + (n-y)*log(1-p)) } # # Estimate some starting values. # m.0 <- x[which.max(y)]; s.0 <- (max(x)-min(x))/4; a.0 <- 0 theta <- c(m.0, s.0, a.0) # # Do the fit. # fit <- nlm(likelihood.log, theta, y, n, x) # # Plot the results. # par(mfrow=c(1,1)) plot(c(min(x),max(x)), c(0,1), main="Data", type="n", xlab="x", ylab="y") curve(gaussian(x, fit$estimate), add=TRUE, col="Red", lwd=2) #$points(x, y/n, pch=19) # # Compute residuals. # mu.hat <- gaussian(x, fit$estimate)
residuals <- (y - n*mu.hat) / sqrt(n * mu.hat * (1-mu.hat))
boxplot(residuals, horizontal=TRUE, xlab="Standardized residuals")
#
# (Compute the variance-covariance matrix in terms of the Hessian
# of the log likelihood at the estimates, etc., etc.)


(Although R has a way of performing these calculations automatically (by means of a custom link function for its glm Generalized Linear Model function), the interface is so cumbersome and so uninstructive I have elected not to use it for this illustration.)

One answer to your question is that least squares is a standard way to fit Gaussians to "some x and y data".

However, your data are special in that they are fractions (probabilities) that must lie in [0,1]. What you are doing is fitting a curve which is for probability densities, so it is unaware of the upper limit. The fitted maximum exceeds 1 and is for a value of x just less than 0. So, from one point of view you fitted an impossible curve. Densities greater than 1 are certainly possible but probabilities greater than 1 certainly are not.

You could scale your probabilities to densities, assuming that each is for a bin width 0.1. However, that would then ignore the fact that your largest value is the largest possible value, so you can't win everything either way.

What you should be doing is a largely open question, as you would need to tell us more about your problem and your data to get better advice. But various comments spring to mind:

1. Even if you have a really good reason to be focusing on Gaussians, I would think of other models too, asymmetric as well as symmetric.

2. The data are presented as if already on a scale with time interval 0.1 s. If your data are finer, you should use exact times. If this is the actual time resolution, then be advised that the data are rather coarse for discriminating between possible models. Naturally, there may be excellent scientific or technical reasons for your time interval.

3. One model popular in ecology for proportional abundances has been called the Gaussian logit, which in your terms is a logit for y fitted as a quadratic, with predictors x and x$^2$. Taking each value to represent 100 repetitions, I fitted this model in Stata:

clear
set obs 7
gen x = (_n - 4)/10
mat y = (10,45,90,100,60,10,5)
gen y = y[1,_n]
gen xsq = x^2
glm y x xsq , link(logit) family(binomial 100)
gen py = y/100
twoway function invlogit(_b[_cons]+ _b[x]*x + _b[xsq]*x^2), ra(x) || scatter py x , scheme(s1color) legend(order(1 "predicted" 2 "observed")) xla(-.3(.1).3) ytitle(proportion)


A reference for this model is Jongman, R.H.G., ter Braak, C.J.F. and van Tongerren, O.F.R. 1995. Data analysis in community and landscape ecology. Cambridge: Cambridge University Press.

The model fit is not especially encouraging. Note that in turn a limitation of, or more neutrally a fact about, this model is that it can never predict a maximum value of exactly 1 (or a minimum value of exactly 0). This model does concur with your least-square result in fitting a maximum for a negative x. However, you (apparently) don't have enough data for much more flexible models to be tried comfortably, unless there are parsimonious physically- or psychometrically-based models tailored to the purpose that can be tried.

• Typically for this kind of problem the data is coarse. Many times a gaussian function is fitted and the mean is called the point of subjective simultaneity. This point often corresponds to an asynchrony in which the flash is presented before the sound and sometimes this is taken as evidence that the perception of the flash needs more time than the sound. So I guess a gaussian is used to easily get some measure of central tendency. Thanks very much for the extensive and clear answer. Sep 24, 2013 at 13:36
• It sounds as if the offset is the parameter of most interest. With a logit model, response the proportion, predictors $x$ and $x^2$, the position of the maximum is $-$coefficient of $x$/(2 * coefficient of $x^2)$. This is just the usual calculus argument that $a + bx + cx^2$ has zero derivative when $x = -b/2c$; we need to check that the value of $x$ really is a maximum. For your data I get $-.03574713$ from a Gaussian logit model. Your R code gives $-0.03636307$. I've given more d.p. than is sensible, but the agreement between slightly different procedures looks good. Sep 24, 2013 at 13:52
• In what's above, "value of $x$" should be "value for $x$ so calculated", i.e. the value of $x$ is the position of the maximum predicted for $y$ (in both the usual and the OP's notation). Sep 24, 2013 at 14:44
• Rounding now more sensibly to 3 d.p. @whuber's method gives $−0.033$ s as compared with $-0.036$ s from the other methods. I would guess that the differences are not worth trying to interpret. Sep 25, 2013 at 12:24
• I agree that the modes are similar in the two fits, but the other two parameters may (ought) to be of interest, too. Their estimates differ distinctly, especially the estimates of the amplitude. A goodness of fit test (or equivalently a comparison of AIC or BIC) will be decisive. This is an important point because it shows that even a small amount of data can discriminate among competing models.
– whuber
Sep 25, 2013 at 14:18