2
$\begingroup$

The problem is from here. Namely, I have set of data points

data = [[90.00, 2.0], [97.40, 5.0], [104.8, 14.0], [112.2, 12.0], [119.6, 11.0], [127.0, 6.0], [134.4, 3.0], [141.8, 1.0], [149.2, 2.0], [156.6, 1.0]]

I have to fit a curve of the form $$\frac{m}{\sigma\sqrt{2\pi}}\cdot e^{-(x-\mu)^2/(2\sigma^2)}$$ to this data.

In the link, Volker Brown got the values

$$m = 405.75796954829985, \mu = 111.86913960269014, \sigma = 11.968861052746961$$

and as I tried my version of Sage to the same to that set of data, I got

$$m = 405.7572223708457, \mu = 111.86914614035226, \sigma = 11.96881745593664$$

What kind of ways there are to measure which one is better fit to that set of data? I heard that sometimes one can use log-likelihoods to select the best model and sometimes Akaike's information criterion is suitable for that. But which one should I choose in this situation? Thusfar, I haven't studied the theory behind the AIC.

Another thing that confuses me is Brown's comment "You need to give some hint for the initial values." Does that mean that there is no standard way to measure how good a particular curve fits to the data?

Or is that just a bug in a previous version of Sage that has been fixed in the version 5.6, 21 January 2013?

$\endgroup$
3
  • 1
    $\begingroup$ How did you measure the values of your data points? What do they represent? What is the purpose of fitting this curve? And why are you concerned about differences in the seventh decimal place in estimates with data that appear to have at most three decimal places of precision and to be limited in quantity? $\endgroup$
    – whuber
    Commented Feb 17, 2013 at 22:15
  • 1
    $\begingroup$ @whuber I don't know where the data is originally from. I found that discussion, tried it and found that my version of Sage outputs different results that the version original poster uses. I'm just curious if there are some methods to evaluate which curve fits better to the data. $\endgroup$ Commented Feb 17, 2013 at 22:21
  • $\begingroup$ These data almost surely are counts associated with numbers representing midpoints of bins. For such data this entire approach would be terribly misguided (and wrong, due to the inclusion of $m$ as a parameter): there are far better ways to fit distributions. Moreover, as a simple plot will show, no Normal distribution will fit these data even approximately well. Your question is a good one but this example is spectacularly bad. $\endgroup$
    – whuber
    Commented Feb 17, 2013 at 23:21

2 Answers 2

4
$\begingroup$

Akaike's information criterion and (log)-likelihood

You are right that the (log-)likelihood and Akaike's information criterion can be used to discriminate between different models.

However, if the two models have the same number of parameters the approaches are equivalent. Akaike's information criterion is defined by (see link) $$ \mathit{AIC} = 2k - 2\ln(L), $$ where $k$ is the number of parameters of the model and $L$ is the likelihood function.

In your case, the number of parameters is $k=2$ for both models such that it is sufficient to consider the likelihood function.

Fitting curves

There are a number of methods to fit curves to data, e.g. least-square fitting or maximum likelihood fitting. Different methods are used in different situations: least-square fitting is appropriate if you want to determine the relationship between two variables. Maximum likelihood fitting is appropriate if you want to determine parameters of a probability distribution.

Because the original problem mentions "I have a problem with fitting my data set with a Gauss (Normal) distribution" I will take a maximum-likelihood approach in the following.

Fitting the parameters (assuming a probability distribution)

The likelihood function is defined as the probability to observe a given sample (your data values) given a set of parameters ($sigma$ and $mu$).

Consider a single data point with value $x$. Assuming that your model is Gaussian, the probability to observe this point is $$p(x|\sigma,\mu)=\frac{1}{\sqrt{2\pi}\sigma} \exp\left(-\frac{1}{2}\frac{(x-\mu)^2}{\sigma^2}\right).$$

Assuming that your samples are independent and are members of a vector $\bf{x}$, the probability to observe them is $$ \begin{align} p({\bf x}|\sigma,\mu)&=\prod_{i=1}^Np(x_i|\sigma,\mu)\\ &=\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{1}{2\sigma^2}\sum_{i=1}^N(x_i-\mu)^2\right). \end{align} $$

We can maximise this probability (the likelihood) by differentiating with respect to $\sigma$ and $\mu$. You can look up the details here.

We thus find $$\hat{\mu} = \overline{x} \equiv \frac{1}{n}\sum_{i=1}^n x_i, \qquad \hat{\sigma}^2 = \frac{1}{n} \sum_{i=1}^n (x_i - \overline{x})^2.$$

I have used the following python code to fit the parameters (assuming that the second value of each tuple is a count).

import numpy as np
data = np.array([[90.00, 2.0], [97.40, 5.0], [104.8, 14.0], [112.2, 12.0], [119.6, 11.0], [127.0, 6.0], [134.4, 3.0], [141.8, 1.0], [149.2, 2.0], [156.6, 1.0]])
#Determine the total number of points
number_of_points = np.sum(data[:,1])
#Determine the mean using a weighted average
mu = np.sum(map(lambda (value, count): value * count, data)) / number_of_points
#Determine the standard deviataion
sigma = np.sqrt(np.sum(map(lambda (value, count): (value - mu)**2 * count, data)) / number_of_points)
print "\\mu = {0}, \\sigma = {1}".format(mu, sigma)

The results are $\mu = 115.056140351, \sigma = 14.1192790625$. Note that we do not need to determine $m$ because it is fixed by the condition that the probability distribution needs to be properly normalised.

Comparing results

The following code evaluates the log-likelihood of the different parameter values

def ll(data, mu, sigma):
    #Determine the total number of points
    number_of_points = np.sum(data[:,1])
    #Determine the variance
    var = np.sum(map(lambda (value, count): (value- mu)**2 * count, data)) / number_of_points

    return - number_of_points*(.5 * np.log(2*np.pi)\
        + np.log(sigma) + var / (2 * sigma ** 2))

print "First result"
print ll(data, 111.86913960269014, 11.968861052746961)
print "Second result"
print ll(data, 111.86914614035226, 11.96881745593664)
print "My result"
print ll(data, mu, sigma)

The results are

First result
-235.552835278
Second result
-235.552923018
My result
-231.789343246

Thus, the result obtained by Volker is marginally better than your SAGE result. However, neither of them maximises the log-likelihood.

Note that all of the above assumes that we are fitting a probability distribution.

$\endgroup$
0
$\begingroup$

You should note that AIC is just a penalized liklihood where the penalty is based on the number of parameters fit. Since the 2 models that you are comparing both fit the same number of parameters (just found different values) they will both receive the same exact penalty. So for your case comparing AIC values will give the same results as comparing liklihoods.

The advantage of AIC is when comparing models with different numbers of parameters (i.e. a linear model with additional covariates, interaction terms, polynomial terms, etc.). The more complicated model will tend to have a higher liklihood, possibly due to random chance. The AIC value looks to see if the increase of liklihood is enough to justify the added complexity. Since your 2 models are equally complex the liklihood part is the only thing that will change.

The comment about the starting values is because many algorythms are sensitive to starting values. You may have found a local max/min or the process may have stopped because the last improvement was small compared to previous steps. Think about the question "Where is your home?", if you are asked this while in your own neighborhood you are likely to give your street name and possibly the exact address, if in a different country you are more likely to only give the city, or possibly even a nearby city that is more likely to be recognized. Similarly the closer your starting values are to the "Truth" the stricter the criteria for stopping.

Have you tried actually plotting your data and the 2 fitted equations? Given that the 2 sets of parameters don't differ until after the 7th significant digit, I doubt that you will see any practical difference. That level of difference is easily explained by differences in starting values rather than bugs or changes in the software (rarely do we ever get the exactly correct answer, just reasonably good approximations).

$\endgroup$
2
  • $\begingroup$ I thought that the curve fitting algorithm runs until it converges to the solution which has at least 17 correct digits. I haven't plotted the equations but I think I won't see difference between curves. $\endgroup$ Commented Feb 18, 2013 at 15:52
  • $\begingroup$ @JaakkoSeppälä, how can the algorithm know if there are 17 correct digits? unless it knows the exact correct answer, but then why iterate at all? These algorithms only look at the degree of improvement in the last step. Since machine precision for a double precision floating point number is about 16 or 17 digits it is not reasonable to expect meaningful steps/comparisons at 17 digits. I am not familiar with Sage, but I would guess that the convergence criteria is probably at around 8 significant digits. $\endgroup$
    – Greg Snow
    Commented Feb 19, 2013 at 2:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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