Data memiliki dua tren; bagaimana cara mengekstrak trendline independen?

34

Saya memiliki satu set data yang tidak dipesan dengan cara tertentu tetapi ketika diplot jelas memiliki dua tren yang berbeda. Regresi linier sederhana tidak akan cukup memadai di sini karena perbedaan yang jelas antara kedua seri. Apakah ada cara sederhana untuk mendapatkan dua trendline linear independen?

Sebagai catatan saya menggunakan Python dan saya cukup nyaman dengan pemrograman dan analisis data, termasuk pembelajaran mesin tetapi bersedia untuk beralih ke R jika benar-benar diperlukan.

masukkan deskripsi gambar di sini

teringat
sumber
6
Jawaban terbaik yang saya miliki sejauh ini adalah untuk mencetak ini di kertas grafik dan menggunakan pensil dan penggaris dan kalkulator ...
jbbiomed
Mungkin Anda bisa menghitung kemiringan berpasangan dan mengelompokkannya menjadi dua "klaster-lereng". Namun ini akan gagal jika Anda memiliki dua tren paralel.
Thomas Jungblut
1
Saya tidak punya pengalaman pribadi dengan itu, tapi saya pikir statsmodels akan layak untuk dicoba. Secara statistik, regresi linier dengan interaksi untuk grup akan memadai (kecuali jika Anda mengatakan Anda memiliki data yang tidak dikelompokkan, dalam hal ini itu sedikit hairier ...)
Matt Parker
1
Sayangnya ini bukan data efek tetapi data penggunaan, dan jelas penggunaan dari dua sistem yang terpisah digabungkan ke dalam set data yang sama. Saya ingin dapat menggambarkan dua pola penggunaan, tetapi saya tidak dapat kembali dan mengingat kembali data karena ini mewakili sekitar 6 tahun nilai informasi yang dikumpulkan oleh klien.
jbbiomed
2
Hanya untuk memastikan: klien Anda tidak memiliki data tambahan yang akan menunjukkan pengukuran mana yang berasal dari populasi yang mana? Ini adalah 100% dari data yang Anda atau klien Anda miliki atau dapat temukan. Juga, 2012 sepertinya pengumpulan data Anda berantakan atau salah satu atau kedua sistem Anda gagal. Membuat saya bertanya-tanya apakah tren naik ke titik itu sangat berarti.
Wayne

Jawaban:

30

Untuk mengatasi masalah Anda, pendekatan yang baik adalah menentukan model probabilistik yang cocok dengan asumsi tentang dataset Anda. Dalam kasus Anda, Anda mungkin menginginkan campuran model regresi linier. Anda dapat membuat model "campuran regressor" yang serupa dengan model campuran gaussian dengan mengaitkan titik data yang berbeda dengan komponen campuran yang berbeda.

Saya telah memasukkan beberapa kode untuk membantu Anda memulai. Kode mengimplementasikan algoritma EM untuk campuran dua regressor (seharusnya relatif mudah diperluas ke campuran yang lebih besar). Kode tampaknya cukup kuat untuk kumpulan data acak. Namun, tidak seperti regresi linier, model campuran memiliki tujuan non-cembung, jadi untuk dataset nyata, Anda mungkin perlu menjalankan beberapa percobaan dengan titik awal acak yang berbeda.

import numpy as np
import matplotlib.pyplot as plt 
import scipy.linalg as lin

#generate some random data
N=100
x=np.random.rand(N,2)
x[:,1]=1

w=np.random.rand(2,2)
y=np.zeros(N)

n=int(np.random.rand()*N)
y[:n]=np.dot(x[:n,:],w[0,:])+np.random.normal(size=n)*.01
y[n:]=np.dot(x[n:,:],w[1,:])+np.random.normal(size=N-n)*.01


rx=np.ones( (100,2) )
r=np.arange(0,1,.01)
rx[:,0]=r

#plot the random dataset
plt.plot(x[:,0],y,'.b')
plt.plot(r,np.dot(rx,w[0,:]),':k',linewidth=2)
plt.plot(r,np.dot(rx,w[1,:]),':k',linewidth=2)

# regularization parameter for the regression weights
lam=.01

def em():
    # mixture weights
    rpi=np.zeros( (2) )+.5

    # expected mixture weights for each data point
    pi=np.zeros( (len(x),2) )+.5

    #the regression weights
    w1=np.random.rand(2)
    w2=np.random.rand(2)

    #precision term for the probability of the data under the regression function 
    eta=100

    for _ in xrange(100):
        if 0:
            plt.plot(r,np.dot(rx,w1),'-r',alpha=.5)
            plt.plot(r,np.dot(rx,w2),'-g',alpha=.5)

        #compute lhood for each data point
        err1=y-np.dot(x,w1)
        err2=y-np.dot(x,w2)
        prbs=np.zeros( (len(y),2) )
        prbs[:,0]=-.5*eta*err1**2
        prbs[:,1]=-.5*eta*err2**2

        #compute expected mixture weights
        pi=np.tile(rpi,(len(x),1))*np.exp(prbs)
        pi/=np.tile(np.sum(pi,1),(2,1)).T

        #max with respect to the mixture probabilities
        rpi=np.sum(pi,0)
        rpi/=np.sum(rpi)

        #max with respect to the regression weights
        pi1x=np.tile(pi[:,0],(2,1)).T*x
        xp1=np.dot(pi1x.T,x)+np.eye(2)*lam/eta
        yp1=np.dot(pi1x.T,y)
        w1=lin.solve(xp1,yp1)

        pi2x=np.tile(pi[:,1],(2,1)).T*x
        xp2=np.dot(pi2x.T,x)+np.eye(2)*lam/eta
        yp2=np.dot(pi[:,1]*y,x)
        w2=lin.solve(xp2,yp2)

        #max wrt the precision term
        eta=np.sum(pi)/np.sum(-prbs/eta*pi)

        #objective function - unstable as the pi's become concentrated on a single component
        obj=np.sum(prbs*pi)-np.sum(pi[pi>1e-50]*np.log(pi[pi>1e-50]))+np.sum(pi*np.log(np.tile(rpi,(len(x),1))))+np.log(eta)*np.sum(pi)
        print obj,eta,rpi,w1,w2

        try:
            if np.isnan(obj): break
            if np.abs(obj-oldobj)<1e-2: break
        except:
            pass

        oldobj=obj

    return w1,w2


#run the em algorithm and plot the solution
rw1,rw2=em()
plt.plot(r,np.dot(rx,rw1),'-r')
plt.plot(r,np.dot(rx,rw2),'-g')

plt.show()
pengguna1149913
sumber
25

Di tempat lain di utas ini, user1149913 memberikan saran yang bagus (mendefinisikan model probabilistik) dan kode untuk pendekatan yang kuat (estimasi EM). Masih ada dua masalah yang harus diatasi:

  1. Bagaimana cara mengatasi penyimpangan dari model probabilistik (yang sangat jelas dalam data 2011-2012 dan agak jelas dalam undulasi dari titik yang kurang miring).

  2. Cara mengidentifikasi nilai awal yang baik untuk algoritma EM (atau algoritma lainnya).

Untuk mengatasi # 2, pertimbangkan untuk menggunakan transformasi Hough . Ini adalah algoritma pendeteksian fitur yang, untuk menemukan peregangan linier fitur, dapat secara efisien dihitung sebagai transformasi Radon .

xyx,ydalam transformasi Hough. Ketika fitur dalam plot asli jatuh di sepanjang garis umum, atau cukup dekat dengan satu, maka koleksi kurva yang mereka hasilkan dalam transformasi Hough cenderung memiliki persimpangan umum yang sesuai dengan garis umum itu. Dengan menemukan titik-titik dengan intensitas terbesar dalam transformasi Hough, kita dapat membacakan solusi yang baik untuk masalah aslinya.

Untuk memulai dengan data ini, saya pertama-tama memangkas hal-hal tambahan (sumbu, tanda centang, dan label) dan untuk ukuran yang baik memangkas titik-titik yang jelas-jelas terpencil di kanan bawah dan menaburkan di sepanjang sumbu bawah. (Ketika hal-hal itu tidak dipotong, prosedur masih bekerja dengan baik, tetapi ia juga mendeteksi sumbu, bingkai, urutan linier kutu, urutan linier label, dan bahkan titik-titik yang terletak secara sporadis pada sumbu bawah!)

img = Import["http://i.stack.imgur.com/SkEm3.png"]
i = ColorNegate[Binarize[img]]
crop2 = ImageCrop[ImageCrop[i, {694, 531}, {Left, Bottom}], {565, 467}, {Right, Top}]

(Ini dan sisa kode ada di Mathematica .)

Cropped image

Untuk setiap titik dalam gambar ini sesuai dengan rentang sempit kurva dalam transformasi Hough, terlihat di sini. Mereka adalah gelombang sinus:

hough2 = Radon[crop2, Method -> "Hough"]  // ImageAdjust

Hough transform

Ini membuat secara visual memanifestasikan arti di mana pertanyaannya adalah masalah pengelompokan garis : transformasi Hough menguranginya menjadi masalah pengelompokan titik , di mana kita dapat menerapkan metode pengelompokan apa pun yang kita suka.

Dalam hal ini, pengelompokan sangat jelas sehingga pasca-pemrosesan sederhana dari transformasi Hough sudah memadai. Untuk mengidentifikasi lokasi dengan intensitas terbesar dalam transformasi, saya meningkatkan kontras dan mengaburkan transformasi pada radius sekitar 1%: itu sebanding dengan diameter titik plot pada gambar asli.

blur = ImageAdjust[Blur[ImageAdjust[hough2, {1, 0}], 8]]

Blurred transform

Thresholding hasilnya mempersempitnya menjadi dua gumpalan kecil yang sentroidnya cukup mengidentifikasi titik-titik intensitas terbesar: ini memperkirakan garis yang pas.

comp = MorphologicalComponents[blur, 0.777]) // Colorize

0,777

Thresholded binarized transform

Sisi kiri gambar sesuai dengan arah 0 derajat (horizontal) dan, seperti yang kita lihat dari kiri ke kanan, sudut itu meningkat secara linear hingga 180 derajat. Interpolasi, saya menghitung bahwa kedua gumpalan itu berpusat pada 19 dan 57,1 derajat, masing-masing. Kita juga dapat membacakan intersepsi dari posisi vertikal gumpalan. Informasi ini menghasilkan kecocokan awal:

width = ImageDimensions[blur][[1]];
slopes =  Module[{x, y, z}, ComponentMeasurements[comp, "Centroid"] /. 
          Rule[x_, {y_, z_}] :>  Round[((y - 1/2)/(width - 1))  180., 0.1]
  ]

{19., 57.1}

Dengan cara yang sama, seseorang dapat menghitung intersep yang sesuai dengan lereng ini, memberikan kesesuaian ini:

Fitted lines

(Garis merah sesuai dengan titik merah muda kecil pada gambar sebelumnya dan garis biru sesuai dengan gumpalan aqua yang lebih besar.)

Sebagian besar, pendekatan ini secara otomatis berurusan dengan masalah pertama: penyimpangan dari linearitas mengotori titik-titik intensitas terbesar, tetapi biasanya tidak banyak menggeser mereka. Titik-titik terpencil yang jujur ​​akan menyumbang kebisingan tingkat rendah di seluruh transformasi Hough, yang akan hilang selama prosedur pasca-pemrosesan.

Pada titik ini orang dapat memberikan estimasi ini sebagai nilai awal untuk algoritma EM atau untuk minimizer kemungkinan (yang, jika diberikan estimasi yang baik, akan bertemu dengan cepat). Namun, yang lebih baik adalah menggunakan estimator regresi yang kuat seperti kuadrat terkecil yang berulang secara berulang . Itu mampu memberikan bobot regresi ke setiap titik. Bobot rendah menunjukkan bahwa suatu titik bukan "milik" suatu garis. Manfaatkan bobot ini, jika diinginkan, untuk menetapkan setiap titik ke garis yang sesuai. Kemudian, setelah mengklasifikasikan poin, Anda dapat menggunakan kuadrat terkecil biasa (atau prosedur regresi lainnya) secara terpisah pada dua kelompok poin.

whuber
sumber
1
Gambar memberi tahu seribu kata dan Anda memiliki 5. Ini adalah karya luar biasa dari grafik singkat yang saya buat hanya untuk keperluan pertanyaan ini! Pujian!
jbbiomed
2
Transformasi Hough banyak digunakan di bidang Computer Vision untuk mengidentifikasi garis lurus dalam gambar. Mengapa itu tidak digunakan dalam statistik juga? ;)
Lucas Reis
Itu pertanyaan yang menarik, @Lucas. Dalam banyak aplikasi statistik, meskipun, ada asimetri antarax dan yvariabel: satu dipandang sebagai diukur secara akurat dan yang lainnya dipandang sebagai subjek variasi acak. Untuk aplikasi seperti itu, metode pemrosesan gambar tidak akan memberikan hasil yang cukup tepat (dan di banyak aplikasi, mereka akan gagal memberikan hasil yang bermanfaat sama sekali). Tapi dari waktu ke waktu mungkin ada beberapa tumpang tindih yang bermanfaat antara kedua bidang: yang saya percaya adalah poin yang ingin Anda sampaikan.
whuber
Iya nih. Bayangkan, misalnya, jumlah outlier yang terlibat dalam membandingkan dua gambar untuk mendeteksi jika mereka berasal dari subjek yang sama. Dan, yang paling penting, bayangkan harus melakukannya secara real time. "Kecepatan" adalah faktor yang sangat penting dalam Visi Komputer, dan tidak begitu penting dalam Statistik.
Lucas Reis
@RoyalTerima kasih telah menunjukkan perlunya perbaikan ke salah satu cuplikan kode. Pada saat saya menemukan perubahan yang Anda sarankan itu telah ditolak (dengan benar, karena itu tidak benar, tetapi tidak masalah bahwa: Saya bersyukur Anda memperhatikan ada kesalahan). Saya memperbaikinya dengan menghapus referensi rotation, yang awalnya diatur ke nol dan karena itu tidak membuat perbedaan.
whuber
15

Saya menemukan pertanyaan ini ditautkan dengan pertanyaan lain . Saya sebenarnya melakukan penelitian akademis tentang masalah seperti ini. Silakan periksa jawaban saya "Least square root" pas? Metode pas dengan beberapa minimum untuk detail lebih lanjut.

Pendekatan berbasis transformasi Hough whuber adalah solusi yang sangat baik untuk skenario sederhana seperti yang Anda berikan. Saya mengerjakan skenario dengan data yang lebih kompleks, seperti ini:

data association problem - candy data set

Rekan penulis dan saya menyatakan ini sebagai masalah "asosiasi data". Ketika Anda mencoba menyelesaikannya, masalah utama biasanya kombinatorial karena jumlah eksponensial dari kombinasi data yang mungkin.

Kami memiliki publikasi " Campuran Tumpang tindih Proses Gaussian untuk masalah asosiasi data " di mana kami mendekati masalah umum kurva N dengan teknik berulang, memberikan hasil yang sangat baik. Anda dapat menemukan kode Matlab tertaut di kertas.

[Perbarui] Implementasi Python dari teknik OMGP dapat ditemukan di perpustakaan GPClust .

Saya memiliki makalah lain di mana kami mengendurkan masalah untuk mendapatkan masalah optimasi cembung, tetapi belum diterima untuk publikasi. Ini khusus untuk 2 kurva, sehingga akan bekerja dengan sempurna pada data Anda. Beri tahu saya jika Anda tertarik.

Steven
sumber
1
Saya sedih melihat bahwa lebih dari dua tahun tidak ada orang lain yang telah mengangkat jawaban asli dan berharga ini. Sementara itu, apakah makalah terakhir yang Anda sebutkan diterima?
whuber
1
Makalah itu memang telah diterima, hanya beberapa bulan yang lalu. Anda dapat mengunduhnya di sini gtas.unican.es/pub/378 . Ini sebenarnya masalah yang sangat jarang (yang mungkin menjelaskan kurangnya popularitas), tetapi kami masih berhasil menemukan beberapa aplikasi yang menarik. Lihatlah eksperimen di akhir makalah jika Anda mau.
Steven
2

user1149913 memiliki jawaban yang sangat baik (+1), tetapi bagi saya tampaknya pengumpulan data Anda berantakan pada akhir 2011, jadi Anda harus memotong bagian data Anda, dan kemudian masih menjalankan beberapa kali dengan acak berbeda mulai koefisien untuk melihat apa yang Anda dapatkan.

Salah satu cara mudah untuk melakukan sesuatu adalah dengan memisahkan data Anda menjadi dua set dengan mata, kemudian menggunakan teknik model linier apa pun yang biasa Anda gunakan. Dalam R, itu akan menjadi lmfungsinya.

Atau paskan dua garis dengan mata. Dalam R Anda akan menggunakannya ablineuntuk melakukan ini.

Data bercampur aduk, memiliki outlier, dan berantakan pada akhirnya, namun dengan mata memiliki dua garis yang cukup jelas, jadi saya tidak yakin metode mewah itu sepadan.

Wayne
sumber