Как я могу использовать метод наименьших квадратов в python, используя данные, которые являются только верхним пределом?

Я пытаюсь выполнить наименьшие квадраты, подходящие в python для известной функции с тремя переменными. Я могу выполнить эту задачу для случайно сгенерированных данных с ошибками, но фактические данные, которые мне нужно подогнать, включают некоторые точки данных, которые являются верхними пределами значений. Функция описывает поток как функцию длины волны, но в некоторых случаях поток, измеренный на заданной длине волны, является не абсолютным значением с ошибкой, а скорее максимальным значением потока, при этом реальное значение может быть любым ниже этого значения вплоть до нуля. .

Есть ли способ сообщить задаче подбора, что некоторые точки данных являются верхними пределами? Кроме того, мне нужно сделать это для ряда наборов данных, и количество точек данных, которые могут быть верхними пределами, различается для каждого из них, поэтому возможность делать это автоматически была бы полезной, но не обязательной.

Я извиняюсь, если что-то из этого неясно, я постараюсь объяснить это более ясно, если это необходимо.

Код, который я использую для соответствия моим данным, приведен ниже.

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

python model-fitting least-squares
person Stargazer_Scot    schedule 08.01.2014    source источник
comment
Имеет ли смысл просто нормально относиться к верхним предельным точкам, но наказывать любое попадание в нефизическую область, например. возвращая в этом случае очень большое значение из residual?   -  person Benjamin Bannier    schedule 08.01.2014
comment
Подгонка не обязательно становится нефизической при подгонке с верхним пределом, но она может ограничивать хвостовую часть наклона (в основном это модифицированная функция черного тела). Без этого ограничения соответствие все еще может быть реальным, но не таким хорошим. Таким образом, обработка верхних предельных точек обычно приводит к плохой подгонке, как в случае с точкой данных с большой ошибкой, а не к предоставлению хотя бы некоторого полезного ввода для функции подбора. Надеюсь, мой ответ имеет смысл.   -  person Stargazer_Scot    schedule 08.01.2014
comment
Я не говорю, что вопрос здесь не по теме, но если вы хотите сначала обсудить методологию того, что вы пытаетесь сделать, stats.stackexchange.com будет хорошим местом.   -  person NPE    schedule 08.01.2014

Ответы (1)


arrow_upward
6
arrow_downward

Scipy.optimize.leastsq — это удобный способ подгонки данных, но основная работа заключается в минимизации функции. Scipy.optimize содержит множество функций минимизации, некоторые из которых способны обрабатывать ограничения. Здесь я объясняю с fmin_slsqp, что я знаю, возможно, другие тоже могут; см. документ Scipy.optimize

fmin_slsqp требуется функция для минимизации и начальное значение параметра. Минимизируемая функция представляет собой сумму квадратов остатков. Для параметров я сначала выполняю традиционную подгонку методом наименьших квадратов и использую результат в качестве начального значения для задачи минимизации с ограничениями. Тогда есть несколько способов наложить ограничения (см. документ); проще параметры f_ieqcons: для этого требуется функция, которая возвращает массив, значения которого всегда должны быть положительными (это ограничения). Здесь функция возвращает положительные значения, если для всех точек максимальных значений функция подгонки находится ниже точки.

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()

В этом примере я подогнал воображаемую выборку точек на параболе. Вот результат без ограничения и с ограничением (красный крест слева): Результаты подгонки

Я надеюсь, что это подойдет для вашего образца данных; в противном случае отправьте один из ваших файлов данных, чтобы мы могли попробовать с реальными данными. Я знаю, что мой пример не учитывает планки ошибок в данных, но вы можете легко справиться с ними, изменив функцию остатков.

person JPG    schedule 09.01.2014
comment
Это сработало отлично. Большое спасибо! Я потратил две недели, пытаясь понять это для себя, и Интернет решает это менее чем за день! Теперь мне просто нужно взять этот код и автоматизировать его для всех моих данных, хотя, надеюсь, это будет самая простая часть! - person Stargazer_Scot; 09.01.2014