2011-10-31 20 views
6

Tôi cần so sánh một số dữ liệu lý thuyết với dữ liệu thực trong python. Dữ liệu lý thuyết đến từ việc giải quyết một phương trình. Để cải thiện so sánh, tôi muốn loại bỏ các điểm dữ liệu nằm xa đường cong lý thuyết. Ý tôi là, tôi muốn loại bỏ các điểm bên dưới và trên các đường đứt nét màu đỏ trong hình (được làm bằng matplotlib). Data points and theoretical curvesXóa các điểm dữ liệu bên dưới một đường cong với python

Cả hai đường cong lý thuyết và điểm dữ liệu đều là các mảng có độ dài khác nhau.

tôi có thể cố gắng để loại bỏ các điểm theo một cách xấp xỉ mắt, ví dụ: các điểm trên đầu tiên có thể được phát hiện bằng:

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')])  

Nhưng tôi muốn sử dụng một ít cách xấp xỉ mắt.

Vì vậy, bất kỳ ai cũng có thể giúp tôi tìm ra cách dễ dàng để xóa các điểm có vấn đề không?

Cảm ơn bạn!

+18

Chỉ cần hoàn toàn từ quan điểm khoa học, tôi sẽ không loại bỏ các điểm trừ khi có lý do hợp lệ R that RÀNG mà bạn cho rằng chúng sai. Bạn có đủ dữ liệu mà các điểm bên ngoài sẽ không có bất kỳ ảnh hưởng nào đến mức phù hợp, vì vậy việc xóa chúng chỉ phục vụ để làm cho biểu đồ trông đẹp, không phục vụ bất kỳ mục đích khoa học nào. – NickLH

+0

Bạn nói đúng, nhưng tôi đã được bảo. –

Trả lời

1

Chỉ cần nhìn vào sự khác biệt giữa đường cong màu đỏ và các điểm, nếu nó lớn hơn sự khác biệt giữa đường cong màu đỏ và đường cong màu đỏ đứt nét, hãy loại bỏ nó.

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

Nhưng hãy bình luận từ NickLH nghiêm túc. Dữ liệu của bạn trông khá tốt mà không có bất kỳ bộ lọc nào, "ngoại lệ" của bạn đều có lỗi rất lớn và sẽ không ảnh hưởng nhiều đến mức phù hợp.

+1

Nó là một ide tốt, nhưng tôi thấy khó để tính toán sự khác biệt giữa đường cong màu đỏ và các điểm, vì cả hai mảng có độ dài khác nhau. Người ta có thể sử dụng nội suy để tạo ra một mảng đường cong màu đỏ với chiều dài mảng điểm. –

+0

Các red_curves có thể được thực hiện với một hàm, chỉ là các giá trị x có liên quan trong đó. – tillsten

4

Đây có thể là quá mức cần thiết và dựa trên nhận xét của bạn

Cả đường cong lý thuyết và các điểm dữ liệu là mảng của độ dài khác nhau.

tôi sẽ làm như sau:

  1. Truncate dữ liệu thiết lập để cho giá trị x của nó nằm trong tối đa và giá trị tối thiểu của tập lý thuyết.
  2. Nội suy đường cong lý thuyết bằng cách sử dụng scipy.interpolate.interp1d và các giá trị x dữ liệu cắt ngắn ở trên. Lý do cho bước (1) là để đáp ứng các ràng buộc của interp1d.
  3. Sử dụng numpy.where để tìm giá trị y dữ liệu nằm ngoài phạm vi giá trị lý thuyết có thể chấp nhận được.
  4. DONT loại bỏ các giá trị này, như đã được đề xuất trong các nhận xét và các câu trả lời khác. Nếu bạn muốn cho rõ ràng, chỉ cho họ ra bằng cách vẽ 'inliners' một màu và 'ngoại lệ' một màu khác.

Đây là tập lệnh gần với những gì bạn đang tìm kiếm, tôi nghĩ vậy. Nó hy vọng sẽ giúp bạn thực hiện những gì bạn muốn:

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

Con số này là kết quả: enter image description here

0

Hoặc bạn có thể sử dụng numpy.where() để xác định cặp xy đáp ứng tiêu chí âm mưu của bạn, hoặc có lẽ liệt kê để làm khá nhiều điều tương tự.Ví dụ:

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 

tôi chắc chắn rằng bạn có thể thay đổi các điều kiện để '2' và '5' trong ví dụ trên là các chức năng của đường cong của bạn

4

Cuối cùng tôi sử dụng một số các Yann mã:

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 

Sử dụng những chức năng cuối cùng tôi đã có thể có được:

Data selection

Cảm ơn tất cả vì sự giúp đỡ của bạn.

+2

+1 Để hiển thị cho chúng tôi giải pháp của bạn và cũng để giữ các điểm xa trên biểu đồ. – NickLH

Các vấn đề liên quan