Удалить точки данных под кривой с помощью Python

Question

Мне нужно сравнить некоторые теоретические данные с реальными данными в python. Теоретические данные получены из решения уравнения. Чтобы улучшить сравнительные данные, я бы хотел удалить точки данных, которые находятся далеко от теоретической кривой. Я имею в виду, я хочу удалить точки ниже и выше красных пунктирных линий на рисунке (сделано с помощью matplotlib).

Как теоретические кривые, так и точки данных являются массивами различной длины.

Я могу попытаться удалить точки примерно так, например: первую верхнюю точку можно определить с помощью:

data2[(data2.redshift<0.4)&data2.dmodulus>1]
rec.array([('1997o', 0.374, 1.0203223485103787, 0.44354759972859786)], dtype=[('SN_name', '|S10'), ('redshift', '<f8'), ('dmodulus', '<f8'), ('dmodulus_error', '<f8')])    

Но я бы хотел использовать менее грубый взгляд.

Итак, кто-нибудь может помочь мне найти простой способ удаления проблемных точек?

Спасибо!

Answer 1

В конце я использую код Yann:

def theoryYatDataX(theoryX,theoryY,dataX):
'''For every dataX point, find interpolated theoryY value. theoryx needed
for interpolation.'''
f = interpolate.interp1d(theoryX,theoryY)
return f(dataX[np.where(dataX<np.max(theoryX))])

def findOutlierSet(data,interpTheoryY,theoryErr):
    '''Find where theoryY-theoryErr < dataY theoryY+theoryErr and return
    valid indicies.'''

    up = np.where(data.dmodulus > (interpTheoryY+theoryErr))
    low = np.where(data.dmodulus < (interpTheoryY-theoryErr))
    # join all the index together in a flat array
    out = np.hstack([up,low]).ravel()

    index = np.array(np.ones(len(data),dtype=bool))
    index[out]=False

    datain = data[index]
    dataout = data[out]

    return datain, dataout

def selectdata(data,theoryX,theoryY):
    """
    Data selection: z<1 and +-0.5 LFLRW separation
    """
    # Select data with redshift z<1
    data1 = data[data.redshift < 1]

    # From modulus to light distance:
    data1.dmodulus, data1.dmodulus_error = modulus2distance(data1.dmodulus,data1.dmodulus_error)

    # redshift data order
    data1.sort(order='redshift')

    # Outliers: distance to LFLRW curve bigger than +-0.5
    theoryErr = 0.5
    # Theory curve Interpolation to get the same points as data
    interpy = theoryYatDataX(theoryX,theoryY,data1.redshift)

    datain, dataout = findOutlierSet(data1,interpy,theoryErr)
    return datain, dataout

Используя эти функции, я наконец могу получить:

Спасибо всем заваша помощь.

Answer 2

Это может быть излишним и основано на вашем комментарии

Как теоретические кривые, так и точки данных являются массивами разной длины.

Я бы сделал следующее:

1. Обрезать набор данных так, чтобы его значения x находились в пределах значений max и min теоретического набора.
2. Интерполируйте теоретическую кривую, используя scipy.interpolate.interp1d и приведенные выше усеченные значения данных x. Причина для шага (1) состоит в том, чтобы удовлетворить ограничения interp1d.
3. Используйте numpy.where, чтобы найти значения y данных, выходящие за пределы допустимых теоретических значений.
4. НЕ отбрасывайте эти значения, как это было предложено в комментариях и других ответах. Если вы хотите для ясности, укажите их, нанеся на график «вкладыши» одного цвета, а «выбросы» - другого цвета.

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

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

# make up data
def makeUpData():
    '''Make many more data points (x,y,yerr) than theory (x,y),
    with theory yerr corresponding to a constant "sigma" in y, 
    about x,y value'''
    NX= 150
    dataX = (np.random.rand(NX)*1.1)**2
    dataY = (1.5*dataX+np.random.rand(NX)**2)*dataX
    dataErr = np.random.rand(NX)*dataX*1.3
    theoryX = np.arange(0,1,0.1)
    theoryY = theoryX*theoryX*1.5
    theoryErr = 0.5
    return dataX,dataY,dataErr,theoryX,theoryY,theoryErr

def makeSameXrange(theoryX,dataX,dataY):
    '''
    Truncate the dataX and dataY ranges so that dataX min and max are with in
    the max and min of theoryX.
    '''
    minT,maxT = theoryX.min(),theoryX.max()
    goodIdxMax = np.where(dataX<maxT)
    goodIdxMin = np.where(dataX[goodIdxMax]>minT)
    return (dataX[goodIdxMax])[goodIdxMin],(dataY[goodIdxMax])[goodIdxMin]

# take 'theory' and get values at every 'data' x point
def theoryYatDataX(theoryX,theoryY,dataX):
    '''For every dataX point, find interpolated thoeryY value. theoryx needed
    for interpolation.'''
    f = interpolate.interp1d(theoryX,theoryY)
    return f(dataX[np.where(dataX<np.max(theoryX))])

# collect valid points
def findInlierSet(dataX,dataY,interpTheoryY,thoeryErr):
    '''Find where theoryY-theoryErr < dataY theoryY+theoryErr and return
    valid indicies.'''
    withinUpper = np.where(dataY<(interpTheoryY+theoryErr))
    withinLower = np.where(dataY[withinUpper]
                    >(interpTheoryY[withinUpper]-theoryErr))
    return (dataX[withinUpper])[withinLower],(dataY[withinUpper])[withinLower]

def findOutlierSet(dataX,dataY,interpTheoryY,thoeryErr):
    '''Find where theoryY-theoryErr < dataY theoryY+theoryErr and return
    valid indicies.'''
    withinUpper = np.where(dataY>(interpTheoryY+theoryErr))
    withinLower = np.where(dataY<(interpTheoryY-theoryErr))
    return (dataX[withinUpper],dataY[withinUpper],
            dataX[withinLower],dataY[withinLower])
if __name__ == "__main__":

    dataX,dataY,dataErr,theoryX,theoryY,theoryErr = makeUpData()

    TruncDataX,TruncDataY = makeSameXrange(theoryX,dataX,dataY)

    interpTheoryY = theoryYatDataX(theoryX,theoryY,TruncDataX)

    inDataX,inDataY = findInlierSet(TruncDataX,TruncDataY,interpTheoryY,
                                    theoryErr)

    outUpX,outUpY,outDownX,outDownY = findOutlierSet(TruncDataX,
                                                    TruncDataY,
                                                    interpTheoryY,
                                                    theoryErr)
    #print inlierIndex
    fig = plt.figure()
    ax = fig.add_subplot(211)

    ax.errorbar(dataX,dataY,dataErr,fmt='.',color='k')
    ax.plot(theoryX,theoryY,'r-')
    ax.plot(theoryX,theoryY+theoryErr,'r--')
    ax.plot(theoryX,theoryY-theoryErr,'r--')
    ax.set_xlim(0,1.4)
    ax.set_ylim(-.5,3)
    ax = fig.add_subplot(212)

    ax.plot(inDataX,inDataY,'ko')
    ax.plot(outUpX,outUpY,'bo')
    ax.plot(outDownX,outDownY,'ro')
    ax.plot(theoryX,theoryY,'r-')
    ax.plot(theoryX,theoryY+theoryErr,'r--')
    ax.plot(theoryX,theoryY-theoryErr,'r--')
    ax.set_xlim(0,1.4)
    ax.set_ylim(-.5,3)
    fig.savefig('findInliers.png')

Эта цифра является результатом:

Answer 3

Просто посмотрите на разницу между красной кривой и точками, если она больше, чем разница между красной кривой и пунктирной красной кривой, удалите ее.

diff=np.abs(points-red_curve)
index= (diff>(dashed_curve-redcurve))
filtered=points[index]

Но, пожалуйста, отнеситесь к комментариям NickLH серьезно. Ваши данные выглядят довольно хорошо без какой-либо фильтрации, все ваши «outlieres» имеют очень большую ошибку и не сильно повлияют на подгонку.

Answer 4

Либо вы можете использовать numpy.where (), чтобы определить, какие пары xy соответствуют вашим критериям построения, либо, возможно, перечислить, чтобы сделать почти то же самое. Пример:

x_list = [ 1,  2,  3,  4,  5,  6 ]
y_list = ['f','o','o','b','a','r']

result = [y_list[i] for i, x in enumerate(x_list) if 2 <= x < 5]

print result

Я уверен, что вы можете изменить условия так, чтобы '2' и '5' в приведенном выше примере были функциями ваших кривых