How can I do a least squares fit in python, using data that is only an upper limit?

Question

I am trying to perform a least squares fit in python to a known function with three variables. I am able to complete this task for randomly generated data with errors, but the actual data that I need to fit includes some data points that are upper limits on the values. The function describes the flux as a function of wavelength, but in some cases the flux measured at the given wavelength is not an absolute value with an error but rather a maximum value for the flux, with the real value being anything below that down to zero.

Is there some way of telling the fitting task that some data points are upper limits? Additionally, I have to do this for a number of data sets, and the number of data points which could be upper limits is different for each one, so being able to do this automatically would be beneficial but not a necessity.

I apologise if any of this is unclear, I will endeavour to explain it more clearly if it is needed.

The code I am using to fit my data is included below.

import numpy as np
from scipy.optimize import leastsq
import math as math
import matplotlib.pyplot as plt


def f_all(x,p):
    return np.exp(p[0])/((x**(3+p[1]))*((np.exp(14404.5/((x*1000000)*p[2])))-1))

def residual(p,y,x,error):
    err=(y-(f_all(x,p)))/error
    return err


p0=[-30,2.0,35.0]

data=np.genfromtxt("./Data_Files/Object_001")
wavelength=data[:,0]
flux=data[:,1]
errors=data[:,2]

p,cov,infodict,mesg,ier=leastsq(residual, p0, args = (flux, wavelength, errors), full_output=True)

print p

Answer 1

Scipy.optimize.leastsq is a convenient way to fit data, but the work underneath is the minimization of a function. Scipy.optimize contains many minimization functions, some of then having the capacity of handling constraints. Here I explain with fmin_slsqp which I know, perhaps the others can do also; see Scipy.optimize doc

fmin_slsqp requires a function to minimize and an initial value for the parameter. The function to minimize is the sum of the square of the residuals. For the parameters, I perform first a traditional leastsq fit and use the result as an initial value for the constrained minimization problem. Then there are several ways to impose constraints (see doc); the simpler is the f_ieqcons parameters: it requires a function which returns an array whose values must always be positive (that's the constraints). Here the function returns positive values if, for all maximal values points, the fit function is below the point.

import numpy
import scipy.optimize as scimin
import matplotlib.pyplot as mpl

datax=numpy.array([1,2,3,4,5]) # data coordinates
datay=numpy.array([2.95,6.03,11.2,17.7,26.8])
constraintmaxx=numpy.array([0]) # list of maximum constraints
constraintmaxy=numpy.array([1.2])

# least square fit without constraints
def fitfunc(x,p): # model $f(x)=a x^2+c
    a,c=p
    return c+a*x**2
def residuals(p): # array of residuals
    return datay-fitfunc(datax,p)
p0=[1,2] # initial parameters guess
pwithout,cov,infodict,mesg,ier=scimin.leastsq(residuals, p0,full_output=True) #traditionnal least squares fit

# least square fir with constraints
def sum_residuals(p): # the function we want to minimize
    return sum(residuals(p)**2)
def constraints(p): # the constraints: all the values of the returned array will be >=0 at the end
    return constraintmaxy-fitfunc(constraintmaxx,p)
pwith=scimin.fmin_slsqp(sum_residuals,pwithout,f_ieqcons=constraints) # minimization with constraint

# plotting
ax=mpl.figure().add_subplot(1,1,1)
ax.plot(datax,datay,ls="",marker="x",color="blue",mew=2.0,label="Datas")
ax.plot(constraintmaxx,constraintmaxy,ls="",marker="x",color="red",mew=2.0,label="Max points")
morex=numpy.linspace(0,6,100)
ax.plot(morex,fitfunc(morex,pwithout),color="blue",label="Fit without constraints")
ax.plot(morex,fitfunc(morex,pwith),color="red",label="Fit with constraints")
ax.legend(loc=2)
mpl.show()

In this example I fit an imaginary sample of points on a parabola. Here is the result, without and with constraint (the red cross on left):

I hope this will do for your data sample; otherwise, please post one of your data files so that we can try with real data. I know my example does not takes care of error bars on data, but you can easily handle them by modifying the residuals function.