8
$\begingroup$

I am quite new to statistics, so please forgive me for using probably the wrong vocabulary.

I have some data that looks (to me) like a gaussian when plotted.

The data is an extract from a jpeg image. It's a vertical line taken from the image, and only the Red data is used (from RGB).

Here is the full data (27 data points):

> r
 [1] 0.003921569 0.031372549 0.023529412 0.015686275 0.003921569 0.027450980
 [7] 0.003921569 0.015686275 0.031372549 0.105882353 0.305882353 0.490196078
[13] 0.560784314 0.615686275 0.592156863 0.505882353 0.364705882 0.227450980
[19] 0.050980392 0.031372549 0.019607843 0.054901961 0.031372549 0.015686275
[25] 0.027450980 0.003921569 0.011764706

> dput(r)
c(0.00392156862745098, 0.0313725490196078, 0.0235294117647059, 
0.0156862745098039, 0.00392156862745098, 0.0274509803921569, 
0.00392156862745098, 0.0156862745098039, 0.0313725490196078, 
0.105882352941176, 0.305882352941176, 0.490196078431373, 0.56078431372549, 
0.615686274509804, 0.592156862745098, 0.505882352941176, 0.364705882352941, 
0.227450980392157, 0.0509803921568627, 0.0313725490196078, 0.0196078431372549, 
0.0549019607843137, 0.0313725490196078, 0.0156862745098039, 0.0274509803921569, 
0.00392156862745098, 0.0117647058823529)
plot(r)

enter image description here

I would like to find a gaussian that is as close as possible to the plot/data.

I tried with normalmixEM from the R package mixtools.

> fit = normalmixEM(r)

but this seems to try to fit to a mix of two gaussian by default.

I tried to specify that there is only one gaussian using the parameter k:

> fit = normalmixEM(r, k = 1)
Error in normalmix.init(x = x, lambda = lambda, mu = mu, s = sigma, k = k,  : 
  arbmean and arbvar cannot both be FALSE

How can I fit the data?

$\endgroup$
5
  • 2
    $\begingroup$ r is clearly not normal. Its distribution is right skewed (lot of small values close to 0, few large values). You will see this by typing "hist(r)". $\endgroup$
    – Michael M
    Jan 22 '14 at 15:27
  • 1
    $\begingroup$ Tip! Use dput(r) to generate a string that is easily copy'n'pasteable. Now we have to enter the data in r manually... $\endgroup$ Jan 22 '14 at 15:30
  • $\begingroup$ @RasmusBååth thanks, I was looking for that command :) I edited the question. $\endgroup$ Jan 22 '14 at 15:37
  • $\begingroup$ You don't plot a sequence of data values to see the distribution. Does $r$ already represent a frequency or probability density of data values? $\endgroup$ Jan 22 '14 at 15:58
  • $\begingroup$ I have answered this question several times in several contexts. An R solution for a discrete variable like your Index appears at stats.stackexchange.com/a/43004/919; an R solution for a continuous variable is at stats.stackexchange.com/questions/70153/…; and an Excel solution is at stats.stackexchange.com/a/11563/919. $\endgroup$
    – whuber
    Jan 22 '14 at 16:45
10
$\begingroup$

There's a difference between fitting a gaussian distribution and fitting a gaussian density curve. What normalmixEM is doing is the former. What you want is (I guess) the latter.

Fitting a distribution is, roughly speaking, what you'd do if you made a histogram of your data, and tried to see what sort of shape it had. What you're doing, instead, is simply plotting a curve. That curve happens to have a hump in the middle, like what you get by plotting a gaussian density function.

To get what you want, you can use something like optim to fit the curve to your data. The following code will use nonlinear least-squares to find the three parameters giving the best-fitting gaussian curve: m is the gaussian mean, s is the standard deviation, and k is an arbitrary scaling parameter (since the gaussian density is constrained to integrate to 1, whereas your data isn't).

x <- seq_along(r)

f <- function(par)
{
    m <- par[1]
    sd <- par[2]
    k <- par[3]
    rhat <- k * exp(-0.5 * ((x - m)/sd)^2)
    sum((r - rhat)^2)
}

optim(c(15, 2, 1), f, method="BFGS", control=list(reltol=1e-9))
$\endgroup$
1
  • $\begingroup$ Played around with this solution but you beat me :) When playing around I noticed that the initial starting values given to optim mattered a lot, so when using this method make sure to check the fit graphically. $\endgroup$ Jan 22 '14 at 16:17
11
$\begingroup$

I propose to use non-linear least squares for this analysis.

# First present the data in a data-frame
tab <- data.frame(x=seq_along(r), r=r)
#Apply function nls
(res <- nls( r ~ k*exp(-1/2*(x-mu)^2/sigma^2), start=c(mu=15,sigma=5,k=1) , data = tab))

And from the output, I was able to obtain the following fitted "Gaussian curve":

v <- summary(res)$parameters[,"Estimate"]
plot(r~x, data=tab)
plot(function(x) v[3]*exp(-1/2*(x-v[1])^2/v[2]^2),col=2,add=T,xlim=range(tab$x) )

enter image description here

The fit is not amazing... Wouldn't a $x \mapsto \sin(x) / x$ function be a better model?

$\endgroup$
3
  • $\begingroup$ Thanks. I get residual sum-of-squares: 0.01997. I think I am getting exactly the same with the solution from Hong Ooi above. Is the algo the same? Also how do I plot the result of nls? $\endgroup$ Jan 22 '14 at 20:30
  • $\begingroup$ Yeah, the algorithms are the same in the sense that if they work (don't get stuck in some local minimum) they give the same answer. Depending on the value given to method= they can be exactly the same. $\endgroup$ Jan 23 '14 at 9:06
  • 1
    $\begingroup$ I added two lines to generate the plot. $\endgroup$ Jan 23 '14 at 9:17

Not the answer you're looking for? Browse other questions tagged or ask your own question.