Orthogonal Distance Regression in Python for Gaussian fitting

Question

I am trying to do a 1D Gaussian fitting using ODR in Python, but keep getting wrong fitting results.

For simplicity, assume that I have a set of 19 data points. These are the data I want to fit:

import numpy as np
import matplotlib.pyplot as plt
from scipy.odr import ODR, Model, Data

x = np.arange(0,19,1)
y = np.array([5.64998480e+09, 3.03653479e+10, 2.18927521e+11, 6.22541771e+11,
       1.24917901e+12, 2.05145638e+12, 2.92904416e+12, 3.74656109e+12,
       4.36310058e+12, 4.66564452e+12, 4.59701326e+12, 4.17028923e+12,
       3.46549578e+12, 2.60950760e+12, 1.74504950e+12, 9.97650569e+11,
       4.49637554e+11, 1.27693929e+11, 6.10512095e+09])

def func(beta,data):
    x = data
    height,center_x,width = beta
    return height*np.exp(-(((center_x-x)/width)**2/2))

data = Data([x],y)
model = Model(func)
odr = ODR(data, model, [1e12,10,2])
res = odr.run()

plt.figure()
plt.plot(x,y)
plt.plot(x,func(res.beta, x),'-o')

This is what I got

What is wrong with my code?

Thanks!

score 0 · Answer 1 · answered Jun 26 '23 at 15:58

0

I think your guess just isn't close enough. I tried beta0=[4e12, 9, 4] and it converged. You can check the results using res.pprint().

import numpy as np
import matplotlib.pyplot as plt
from scipy.odr import ODR, Model, Data

plt.close("all")

x = np.arange(0,19,1)
y = np.array([5.64998480e+09, 3.03653479e+10, 2.18927521e+11, 6.22541771e+11,
       1.24917901e+12, 2.05145638e+12, 2.92904416e+12, 3.74656109e+12,
       4.36310058e+12, 4.66564452e+12, 4.59701326e+12, 4.17028923e+12,
       3.46549578e+12, 2.60950760e+12, 1.74504950e+12, 9.97650569e+11,
       4.49637554e+11, 1.27693929e+11, 6.10512095e+09])

def func(beta, x):
    height, center_x, width = beta
    return height*np.exp(-(((center_x-x)/width)**2/2))

data = Data(x, y)
model = Model(func)
odr = ODR(data, model, beta0=[4e12, 9, 4])
res = odr.run()
res.pprint()

plt.figure()
plt.plot(x, y)
plt.plot(x, func(res.beta, x),'-o')

answered Jun 26 '23 at 15:58

jared

4,165
1
8
31

When I used scipy.optimize.curve_fit even if I give it initial value [1e8,1,1] it still gives pretty good fitting results almost immediately. Is ODR not suitable for Gaussian fitting? – Chris Lin Jun 26 '23 at 16:07
@ChrisLin This was my first time using `ODR`, so I am not sure. I have used `curve_fit` many times and generally prefer it. – jared Jun 26 '23 at 16:14
[This](https://stackoverflow.com/a/27299335/12131013) answer seems to show that the initial guess is very important for ODR. [This](https://stackoverflow.com/a/62830653/12131013) answer also notes some issues with ODR. – jared Jun 26 '23 at 16:17

score 0 · Answer 2 · answered Jun 27 '23 at 08:34

Your software uses an iterative method of regression starting from "guessed" initial values of the parameters. As judiciously pointed out by jared and others it is important to set the initial values not too far from the unknown correct values.

If you don't know how to set good initial values the method shown below is not iterative and doesn't requires guessed initial values. This method can be used to compute approximate values of the parameters. If the results is not sufficiently accurate one can use theses values as initial values for your software.

Results :

About the method of fitting without initial guess : https://fr.scribd.com/doc/14674814/Regressions-et-equations-integrales

You basically have a repeated answer to a very frequent question on StackOverflow where your paper applies. That's not necessarily a bad thing, but I recommend two (non-exclusive) approaches to these questions. First, you should write a [community wiki](https://stackoverflow.com/help/privileges/community-wiki) with more detail on this class of question. Second, many of these questions should just be closed as duplicates and link to one representative question. — Reinderien, Jun 27 '23 at 14:54

Orthogonal Distance Regression in Python for Gaussian fitting

2 Answers2