1

I'm using scipy.optimize.curve_fit to approximate peaks in my data with Gaussian functions. This works well for strong peaks, but it is more difficult with weaker peaks. However, I think fixing a parameter (say, width of the Gaussian) would help with this. I know I can set initial "estimates" but is there a way that I can easily define a single parameter without changing the function I'm fitting to?

finefoot
  • 9,914
  • 7
  • 59
  • 102
Corey Husic
  • 31
  • 1
  • 3
  • 2
    What have you tried so far? Perhaps you could add an example of your code (see [here](http://stackoverflow.com/help/mcve) for help)? – tmdavison Jul 29 '15 at 16:20
  • 1
    Please add some data and some code, then it will be easier to help you. – Cleb Jul 29 '15 at 17:09
  • 1
    Edit the Gaussian function you're fitting such that the parameter you want to hold constant no longer depends on the input arguments. – ali_m Jul 30 '15 at 01:04
  • 2
    Perhaps [this](http://stackoverflow.com/questions/12208634/fitting-only-one-paramter-of-a-function-with-many-parameters-in-python) helps – Peter9192 May 31 '16 at 08:28

2 Answers2

1

If you want to "fix" a parameter of your fit function, you can just define a new fit function which makes use of the original fit function, yet setting one argument to a fixed value:

custom_gaussian = lambda x, mu: gaussian(x, mu, 0.05)

Here's a complete example of fixing sigma of a Gaussian function to 0.05 (instead of optimal value 0.1). Of course, this doesn't really make sense here because the algorithm has no problem in finding optimal values. Yet, you can see how mu is still found despite the fixed sigma.

import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize 

def gaussian(x, mu, sigma):
    return 1 / sigma / np.sqrt(2 * np.pi) * np.exp(-(x - mu)**2 / 2 / sigma**2)

# Create sample data
x = np.linspace(0, 2, 200)
y = gaussian(x, 1, 0.1) + np.random.rand(*x.shape) - 0.5
plt.plot(x, y, label="sample data")

# Fit with original fit function
popt, _ = scipy.optimize.curve_fit(gaussian, x, y)
plt.plot(x, gaussian(x, *popt), label="gaussian")

# Fit with custom fit function with fixed `sigma`
custom_gaussian = lambda x, mu: gaussian(x, mu, 0.05)
popt, _ = scipy.optimize.curve_fit(custom_gaussian, x, y)
plt.plot(x, custom_gaussian(x, *popt), label="custom_gaussian")

plt.legend()
plt.show()

figure

finefoot
  • 9,914
  • 7
  • 59
  • 102
-1

Hopefully this is helpful. Had to use hax. Curve_fit is pretty strict about what it takes.

import numpy as np
from numpy import random
import scipy as sp
from scipy.optimize import curve_fit
import matplotlib.pyplot as pl

def exp1(t,a1,tau1):
    #A1*exp(-t/t1)
    val=0.
    val=(a1*np.exp(-t/tau1))*np.heaviside(t,0)
    return val

def wrapper(t,*args):

    global hold
    global p0
    wrapperName='exp1(t,'
    for i in range(0, len(hold)):
        if hold[i]:
            wrapperName+=str(p0[i])
        else:
            if i%2==0:
                wrapperName+='args['+str(i)+']'
            else:
                wrapperName+='args'+str(i)+']'
        if i<len(hold):
            wrapperName+=','
    wrapperName+=')'

    return eval(wrapperName)

p0=np.array([1.5,500.])
hold=np.array([0,1])
p1=np.delete(p0,1)

timepoints = np.arange(0.,2000.,20.)
y=exp1(timepoints,1,1000)+np.random.normal(0, .1, size=len(timepoints))

popt, pcov = curve_fit(exp1, timepoints, y, p0=p0)
print 'unheld parameters:', popt, pcov

popt, pcov = curve_fit(wrapper, timepoints, y, p0=p1)
for i in range(0, len(hold)):
    if hold[i]:
        popt=np.insert(popt,i,p0[i])
yfit=exp1(timepoints,popt[0],popt[1])
pl.plot(timepoints,y,timepoints,yfit)
pl.show()
print 'hold parameters:', popt, pcov
Scott Ellis
  • 1
  • 1