Bab II Teori Dasar
II.1 Estimasi Spasial Data spasial adalah data yang memuat informasi lokasi. Misalkan z (si ) , i = 1, 2,…,n, s i ∈ D, adalah data observasi peubah acak Z di lokasi atau koordinat yang dinyatakan dengan vektor s i . Vektor koordinat s i adalah vektor yang menyatakan lokasi pengukuran di suatu lapangan spasial D ⊆ R n . Antara satu lokasi dengan lokasi lainnya terdapat kebergantungan, sehingga data spasial merupakan data dependen karena data spasial dikumpulkan dari lokasi spasial yang berbeda yang mengindikasikan kebergantungan antara pengukuran data dengan lokasi. Koleksi peubah acak {Z(s), s ∈ D} disebut proses spasial. Untuk selanjutnya, peubah acak Z(s) yang berkaitan dengan lokasi disebut peubah acak regional. Suatu proses spasial dikatakan stasioner kuat jika untuk sebarang koordinat titik s1 , s 2 ,...., s n dan untuk sebarang vektor h , yang berdimensi sama dengan s i , berlaku Fs1 ,s2 ,....,sn ( z1 , z2 ,...., zn ) = Fs1 +h ,s2 + h ,....,sn +h ( z1 , z2 ,...., zn )
(2.1)
dengan F adalah fungsi distribusi gabungan yang didefinisikan sebagai Fs1 ,s2 ,....,sn ( z1 , z2 ,...., zn ) = P ( Z (s1 ) ≤ z1 , Z (s 2 ) ≤ z2 ,...., Z (s n ) ≤ zn )
(2.2)
Suatu proses spasial dikatakan mengikuti kestasioneran orde 2 jika memiliki sifatsifat sebagai berikut 1. E[ Z (s)] = m , untuk setiap s ∈ D. Artinya mean dari Z ada dan tidak bergantung pada lokasi . 2. Cov[ Z (s), Z (s + h)] = C(h ) , untuk setiap s, s + h ∈ D. Artinya kovariansi antara Z (s) dan Z (s + h ) di lokasi s dan s + h hanya bergantung pada panjang dan
arah vektor h dan tidak bergantung pada koordinat s.
3
Suatu proses spasial dikatakan bersifat intrinsik jika untuk setiap vektor h dan setiap pasangan s dan s + h di D, berlaku 1. E[ Z (s) − Z (s + h )] = m (h ) = 0 , dan 2.
var[ Z (s) − Z (s + h )] = E [ Z (s) − Z (s + h) ] = 2γ (h) , 2
variansi
dari
selisih
Z (s) dengan Z (s + h ) hanya bergantung pada panjang dan arah vektor h dan tidak
bergantung pada lokasi s. Jika {Z(s), s ∈ D} stasioner kuat atau mengikuti stasioner orde 2 maka ia intrinsik. Namun tidak demikian sebaliknya. Besaran 2 γ (h)
di atas disebut juga variogram, sedangkan γ (h) disebut
semivariogram. Semivariogram adalah perangkat yang digunakan untuk menggambarkan, memodelkan, dan menghitung korelasi spasial antara peubah acak regional Z (s) dan Z (s + h ) .
Semivariogram dapat ditaksir dengan
menggunakan data pengamatan. Penaksir bagi semivariogram disebut sebagai semivariogram eksperimental yang didefinisikan sebagai
1 γˆ (h) = 2 N (h) dengan
N (h ) menyatakan
N (h )
∑ ( z(s i =1
himpunan
i
+ h) − z (si ) )
pasangan
2
koordinat
(2.3)
(si , si + h )
yang
dipisahkan oleh vektor h tanpa memperhatikan urutan, dan N (h) menyatakan banyaknya anggota N (h ) .
Hal-hal
yang
perlu
diperhatikan
dalam
penghitungan
semivariogram
eksperimental adalah sebagai berikut, 1. bila sampel hilang (missing value) dari pola regular, nilai sampel yang hilang tersebut tidak perlu diinterpolasi dengan mengambil meannya atau menggantinya dengan nilai 0, 2. bila data iregular, semivariogram dihitung untuk kelas jarak dengan toleransi tertentu, 3. untuk menghitung semivariogram eksperimental perlu diperhatikan arah dan panjang jarak antara titik sampel, dengan kata lain perlu diperhatikan arah dan panjang vektor h.
4
Jika {Z(s), s ∈ D} stasioner kuat atau mengikuti stationer orde 2, maka berlaku
γ (h ) = C(0) − C(h )
(2.4)
dengan C(h ) adalah kovariansi antara Z dari dua lokasi yang dipisahkan oleh vektor h. Besaran C(h ) disebut juga kovariogram. Sedangkan C(0) adalah suatu konstanta yang merupakan kovariogram untuk h = 0. Kovariogram mempunyai perilaku yang berkebalikan dengan semivariogram. Apabila semivariogram naik untuk suatu jarak pisah tertentu maka kovariogram akan turun pada untuk jarak pisah tersebut. Perhatikan bahwa C(0) merupakan variansi dari peubah region Z(s). Jadi pada h = 0, nilai kovariogram sama dengan variansi. Dari kovariogram C(h) dapat dibentuk struktur korelasi spasial ρ (h ) (atau disebut juga korelogram) dengan rumusan
ρ (h) =
C(h) C(0)
(2.5)
Hubungan antara kovariogram dan semivariogram dapat diilustrasikan dengan gambar berikut
σ
γ ( h)
C (h) h
Gambar II.1. Kurva semivariogram dan kovariogram
Kurva yang turun menggambarkan kovariogram, sedangkan kurva yang naik menggambarkan semivariogram. Terlihat bahwa jika semivariogram naik, maka kovariogram turun. Untuk proses spasial yang stasioner atau mengikuti stasioner orde dua, semivariogramnya terbatas oleh nilai C(0) (variansi), artinya nilai semivariogramnya naik seiring dengan makin besarnya jarak h , dan semakin
5
mendekati nilai C(0) tetapi tidak akan melebihi nilai tersebut. Dengan kata lain lim γ (h) = C(0) . Sedangkan kovariogramnya semakin turun mendekati nol
h →∞
seiring dengan makin besarnya jarak h . Dengan demikian, korelasi antara dua lokasi pengukuran akan semakin berkurang seiring dengan makin besarnya jarak pisah antara kedua lokasi. Menurut Matheron (1972), untuk proses spasial yang stasioner dan intrinsik berlaku lim
γ (h )
h →∞
h
2
=0
(2.6).
Plot semivariogram γˆ (h) terhadap jarak h memberikan plot semivariogram eksperimental. Semivariogram eksperimental dari data biasanya bentuknya tidak beraturan sehingga sulit ditafsirkan dan tidak dapat langsung digunakan. Plot semivariogram eksperimental harus dibuat untuk beberapa arah mata angin yang berbeda, sekurang kurangnya empat arah mata angin yaitu Barat-Timur, UtaraSelatan, Barat Laut – Tenggara, dan Timur Laut – Barat Daya. Apabila plot semivariogram eksperimental untuk keempat arah mata angin tersebut tidak terlalu berbeda secara signifikan maka semivariogram tersebut dikatakan isotropik, artinya nilai semivariogram hanya bergantung pada h
dan tidak
bergantung arah. Selanjutnya nilai semivariogram eksperimental akan dicocokkan dengan model semivariogram teoritis untuk digunakan dalam penaksiran. Beberapa model semivariogram teoritis baku yang sering digunakan diantaranya adalah: 1. model sperikal, 2. model eksponensial, 3. model Gaussian, 4. model linear.
II.1.1 Semivariogram Eksperimental Robust Persamaan 2.3 adalah persamaan semivariogram eksperimental yang merupakan penaksir bagi nilai semivariogram untuk suatu vektor h tertentu. Penaksir semivariogram yang dinyatakan oleh persamaan tersebut tidak robust, atau
6
dengan kata lain sangat mudah terpengaruh oleh data pencilan (outliers) karena menggunakan rata-rata sampel sebagai penaksir bagi ekspektasi. Rata-rata sampel adalah penaksir parameter lokasi yang peka terhadap data pencilan.
Cressie dan Hawkins (1980) mengusulkan suatu rumusan untuk menaksir variogram robust untuk data yang berdistribusi normal. Salah satu penaksir parameter lokasi yang robust adalah median. Maka dengan menggunakan median, diperoleh penaksir robust bagi variogram yaitu 1 ⎡ ⎧ ⎫⎤ 2γˆ (h) = ⎢ median ⎨ Z (si ) − Z (s j ) 2 : (si , s j ) ∈ N (h) ⎬⎥ ⎩ ⎭⎦ ⎣
4
B(h)
(2.7)
dengan N (h) menyatakan himpunan pasangan koordinat (si , s j ) yang dipisahkan
oleh vektor h tanpa memperhatikan urutan, dan B (h) merupakan koreksi bias
( B (h) ≈ 0.457 ).
II.1.2 Kriging Kriging adalah metode estimasi yang memberikan taksiran linear tak bias terbaik (Best Linear Unbiased Estimation) untuk suatu titik atau blok. Terbaik di sini artinya memiliki variansi minimum. Ketepatan estimasi kriging sangat bergantung dari model semivariogram yang dipilih.
Misalkan Z (si ) adalah peubah acak region di titik koordinat si , i = 1, 2, ..., n. Taksiran Zˆ (V ) di suatu daerah V (V dapat berupa sebuah titik atau blok) dinyatakan dengan kombinasi linear dari Z (si ) dengan bobot λ i , yaitu n
Zˆ (V ) = ∑ λi Z (si )
(2.8)
i =1
Untuk kasus khusus di mana V={ s 0 }, maka Zˆ (V ) ditulis sebagai Zˆ (s 0 ) Selisih nilai Z (s 0 ) yang sebenarnya dengan nilai taksiran Zˆ (s 0 ) disebut galat kriging ε (s 0 ) . Galat kriging dinyatakan oleh persamaan
ε (s 0 ) = Zˆ (s 0 ) − Z (s 0 )
7
(2.9)
Penaksir Z (s 0 ) yang diperoleh melaui metode kriging memiliki sifat tak bias (unbiased estimator), dengan kata lain E ( ε (s 0 ) ) = 0 , dan variansi galatnya, Var ( ε (s 0 ) ) , minimum.
II.1.2.1 Ordinary Kriging Proses spasial {Z(s), s ∈ D} diasumsikan stasioner dengan rata-rata µ . Rata-rata tersebut sama dengan µ untuk setiap titik dan juga untuk setiap blok. Dengan kata lain E ( Z (s) ) = µ = E ( Z (s 0 ))
(2.10)
untuk setiap titik s dan s 0 ∈ D.
Dalam ordinary kriging, penaksir bagi Z (s 0 ) dimodelkan sebagai n
Zˆ (s 0 ) = ∑ λi Z (si )
(2.11)
i =1
Ekspektasi dari galat estimasi adalah ⎛ n ⎞ n ⎛ n ⎞ E Zˆ (s 0 ) − Z (s 0 ) = E ⎜ ∑ λi Z (si ) − Z (s 0 ) ⎟ = ∑ λi µ − µ = µ ⎜ ∑ λi − 1⎟ (2.12) ⎝ i =1 ⎠ i =1 ⎝ i =1 ⎠
(
)
(
Agar Zˆ (s 0 ) tak bias, dengan kata lain agar E Zˆ (s 0 ) − Z (s 0 )
)
= 0, maka
n
∑λ i =1
i
haruslah sama dengan 1.
Variansi dari galat estimasi adalah ⎛ n ⎞ Var Zˆ (s 0 ) − Z (s 0 ) = E ⎜ ∑ λi Z (si ) −Z (s 0 ) ⎟ ⎝ i =1 ⎠
(
)
2
(2.13)
Untuk memperoleh taksiran dengan variansi galat yang minimum, maka variansi galat harus diminimumkan. Dengan kata lain, harus dicari koefisien-koefisien 2
λ1 , λ2 ,...., λn sedemikian sehingga fungsi E ⎜⎛ ∑ λi Z (si ) −Z (s 0 ) ⎟⎞ − 2m ⎜⎛ ∑ λi − 1⎟⎞ n
⎝ i =1
8
⎠
n
⎝ i =1
⎠
minimum (parameter m adalah faktor pengali Lagrange untuk memastikan bahwa n
∑λ i =1
i
= 1 ). n
∑λ
Perhatikan bahwa syarat
i =1
= 1 menyebabkan
i
2
2
n ⎛ n ⎞ ⎛ n ⎞ 2 λ Z ( s ) Z ( s ) λ Z ( s ) 2 Z ( s ) − = − 0 ⎟ 0 ∑ λi Z (s i ) + Z (s 0 ) ⎜∑ i i ⎜∑ i i ⎟ i =1 ⎝ i =1 ⎠ ⎝ i =1 ⎠ n
n
n
= ∑∑ λi λ j Z (si ) Z (s j ) − 2Z (s 0 )∑ λi Z (si ) + Z (s 0 ) 2 i =1 j =1
= Z (s 0 ) n
i =1
n
2
n
n
∑ λ − 2Z (s )∑ λ Z (s ) + ∑ λ Z (s ) i =1
i
0
i =1
n
i
i
i =1
i
2
i
n
+ ∑∑ λi λ j Z (si ) Z (s j ) − ∑ λi Z (si ) 2 i =1 j =1
i =1
= 2∑ λi ( Z (s 0 ) 2 − 2Z (s 0 ) Z (si ) + Z (si ) 2 ) 2 n
i =1
− ∑∑ λi λ j ( Z (si ) 2 − 2Z (si ) Z (s j ) + Z (s j ) 2 ) 2 n
n
i =1 j =1
= 2∑ λi ( Z (s 0 ) − Z (si ) ) 2 − ∑∑ λi λ j ( Z (si ) − Z (s j ) ) 2 n
2
i =1
n
n
2
i =1 j =1
Dengan demikian 2
⎛ n ⎞ ⎛ n ⎞ E ⎜ ∑ λi Z (si ) −Z (s 0 ) ⎟ − 2m ⎜ ∑ λi − 1⎟ ⎝ i =1 ⎠ ⎝ i =1 ⎠ n n n ⎛ ⎞ 2 ⎛ n ⎞ 2 = E ⎜ 2∑ λi ( Z (s 0 ) − Z (si ) ) 2 − ∑∑ λi λ j ( Z (si ) − Z (s j ) ) 2 ⎟ − 2m ⎜ ∑ λi − 1⎟ i =1 j =1 ⎝ i =1 ⎠ ⎝ i =1 ⎠ ⎛ ( Z (s ) − Z (s ) )2 ⎞ n ⎛ ( Z (s 0 ) − Z (si ) )2 ⎞ n n i j ⎟ − 2m ⎛ ∑ λi − 1⎞ = 2∑ λi E ⎜ ⎟ − ∑∑ λi λ j E ⎜ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ 2 2 i =1 ⎝ i =1 ⎠ ⎝ ⎠ i =1 j =1 ⎝ ⎠ n
n n n ⎛ n ⎞ = 2∑ λiγ ( s 0 − si ) − ∑∑ λi λ jγ ( si − s j ) − 2m ⎜ ∑ λi − 1⎟ i =1 i =1 j =1 ⎝ i =1 ⎠
(2.14)
Sekarang permasalahannya menjadi menentukan koefisien-koefisien λ1 , λ2 ,...., λn yang meminimumkan bentuk (2.14). Dengan menurunkan bentuk (2.14) terhadap
λ1 , λ2 ,...., λn , dan m, dan menyamakan hasilnya dengan nol, diperoleh sistem persamaan
9
∑ λ γ (s n
i
i =1 n
∑λ
i
i =1
i
− s j ) + m = γ ( s 0 − si ) , i = 1, 2, ...., n
(2.15) =1 n
Jadi koefisien-koefisien λ1 , λ2 ,...., λn untuk penaksir kriging Zˆ (s 0 ) = ∑ λi Z (si ) i =1
dapat
dicari
dengan
menyelesaikan
sistem
persamaan
(2.15)
terhadap
λ1 , λ2 ,...., λn . Variansi dari galat estimasi, atau disebut juga variansi kriging, dapat dinyatakan dengan ⎛ n ⎞ σ ( s0 ) = E ⎜ ∑ λi Z (si ) −Z (s0 ) ⎟ ⎝ i =1 ⎠
2
2 K
n
= ∑ λ i γ ( s 0 − si ) + m i =1
= 2∑ λ i γ ( s 0 − si ) − ∑∑ λ i λ j γ ( si − s j ) n
n
i =1
n
(2.16)
i =1 j=1
Persamaan (2.15) dapat dinyatakan dalam bentuk matriks, yaitu ⎛ γ11 γ12 ⎜ ⎜ γ 21 γ 22 ⎜ ⎜ ⎜ γ n1 γ n 2 ⎜ 1 ⎝1
… γ1n 1 ⎞ ⎛ λ1 ⎞ ⎛ γ 01 ⎞ ⎟⎜ ⎟ ⎜ ⎟ … γ 2n 1 ⎟ ⎜ λ 2 ⎟ ⎜ γ 02 ⎟ ⎟⎜ ⎟ = ⎜ ⎟ ⎟⎜ ⎟ ⎜ ⎟ … γ nn 1 ⎟ ⎜ λ n ⎟ ⎜ γ 0n ⎟ … 1 0 ⎠⎟ ⎝⎜ m ⎠⎟ ⎝⎜ 1 ⎠⎟
(2.17)
dengan γ ij = γ ( si − s j ) , untuk i,j = 0, 1, ..., n.
II.1.2.2 Kriging Mean Tujuan ordinary kriging adalah untuk menaksir nilai variabel regional dengan suatu fungsi linear. Sekarang akan diturunkan fungsi linear untuk menaksir nilai mean yang tidak diketahui pada suatu lapangan spasial. Taksiran bagi mean tersebut dapat dituliskan sebagai berikut n
µˆ = ∑ λ µi Z(si ) i =1
10
(2.18)
Seperti sebelumnya, penaksir tersebut haruslah tak bias dan variansi galatnya minimum. Untuk mendapatkan penaksir yang tak bias, ekspektasi dari galat haruslah sama dengan nol. Jadi haruslah
⎛ n ⎞ E ( µˆ − µ ) = E ⎜ ∑ λ µi Z(si ) − µ ⎟ = 0 ⎝ i =1 ⎠
(2.19)
Karena mean Z(s) adalah µ , maka haruslah n
∑λ i =1
µi
=1
(2.20)
Variansi galat estimasi adalah
⎛ n ⎞ n n Var ( µˆ − µ ) = Var ⎜ ∑ λ µi Z(si ) − µ ⎟ = ∑∑ λ µi λ µjC(si − s j ) ⎝ i =1 ⎠ i =1 j=1
(2.21)
Seperti dalam ordinary kriging, akan dicari bobot α i , i = 1, 2,...., n , yang n
meminimumkan variansi galat estimasi dengan batasan
∑λ i =1
µi
= 1 . Maka
didapatkan bobot λ µi , i = 1, 2,...., n , yang memenuhi merupakan solusi dari sistem persamaan kriging n
∑λ j=1
µj
(2.22)
n
∑λ i =1
C(si − s j ) = mµ , i = 1, 2,...., n
µi
=1
dengan mµ adalah faktor pengali Lagrange. Variansi galat estimasi, atau disebut juga variansi kriging, untuk penaksir ini adalah σ 2K = mµ
(2.23)
II.1.2.3 Sequential Kriging Seperti halnya ordinary kriging, sequential kriging adalah estimator yang meminimumkan variansi galat penaksiran. Kumpulan data (data set) dibagi ke dalam beberapa subset dan tiap subset memungkinkan datum tunggal. Dalam tersedianya suatu data tambahan, sequential estimator memperbaiki estimasi
11
sebelumnya dengan menggunakan bobot linier dari data yang baru dan estimasi sebelumnya di suatu lokasi.
Pendekatan sequential kriging dapat menghilangkan kesukaran numerik seperti pada ordinary kriging untuk menyelesaikan persamaan 2.17. Dengan demikian, kompleksitas komputasi dapat dikurangi. Pada pendekatan ini diperbolehkan menggabungkan suatu kelompok data setiap saat atau secara sekuensial selama estimasi.
Misal Zˆ (s 0 ) ädalah estimasi nilai Z di posisi s0. Berdasarkan data sampel Z ( s1 ) , Z ( s 2 ) ,..., Z ( s n ) di lokasi s1 , s 2 ,..., s n , maka estimasi kriging sequential
dalam k langkah diberikan oleh sistem persamaan : p
Zˆ (1) ( s 0 ) = Zˆ (1) (s) = ∑ θi (1) Z ( si ) i =1
Zˆ (2) ( s 0 ) = Zˆ (2) (s) +
q
∑θ
i = p +1
(2) i
( Z (s ) − Zˆ (s ) ) (1)
i
i
(2.24) Zˆ ( k −1) ( s 0 ) = Zˆ ( k − 2) (s) + Zˆ ( k ) ( s 0 ) = Zˆ ( k −1) (s) +
l
∑θ
i = r +1
( k −1) i
n
∑θ
i = l +1
(k ) i
( Z (s ) − Zˆ (s ) ) ( k − 2)
i
i
( Z (s ) − Zˆ (s ) ) ( k −1)
i
i
Superskrip merepresentasikan langkah atau indeks subsets, Zˆ ( k −1) ( s 0 ) adalah nilai estimasi Z (s 0 ) dengan menggunakan data set ke-(k-1), ( Z (si ) − Zˆi (si )( k −1) ) adalah selisih antara nilai data ke-i dengan nilai estimasi data ke-i berdasarkan (k-1) data set, dan θi ( k ) adalah bobot sequential kriging untuk data set ke-k.
Dengan menggunakan sequential kriging, estimasi pertama Z (s 0 ), Z (s 1 ) dan Z (s 2 ) hanya dengan menggunakan nilai Z (s 1 ) adalah : Zˆ (s 0 ) = λ01Z (s 1 ) , Zˆ (s 1 ) = λ11Z (s1 ) , Zˆ (s 2 ) = λ21Z (s1 )
12
(2.25)
dengan
λ01 =
C (s 0 − s1 ) C (s1 − s1 ) C (s 2 − s1 ) = ρ01 , λ11 = = ρ11 = 1 , λ21 = = ρ 21 C (s1 − s1 ) C (s1 − s1 ) C (s1 − s1 )
(2.26)
Selanjutnya Z (s 0 ) akan diestimasi dengan data tambahan Z (s 2 ) (2) (1) (1) Zˆ ( s 0 ) = Zˆ ( s 0 ) + θ 2 ( Z ( s 2 ) − Zˆ ( s 2 ) )
(2.27)
Bobot ( θ 2 ) diperoleh dari persamaan :
θ2 =
ρ02 − ρ 21 ρ01 1 − ρ 212
(2.28)
Estimasi di persamaan 2.27 kemudian menjadi : (2) (1) (1) Zˆ ( s 0 ) = Zˆ ( s 0 ) + θ 2 ( Z ( s 2 ) − Zˆ ( s 2 ) )
= ( ρ 01 − θ 2 ρ 21 ) Z (s1 ) + θ 2 Z (s 2 )
=
( ρ01 − ρ02 ρ 21 ) Z (s ) + ( ρ02 − ρ01ρ21 ) Z (s 1 − ρ 212
1
1 − ρ 212
Misal diperoleh data tambahan baru Z (s3 ) di lokasi s3 .
2
)
(2.29)
Langkah kriging
sekuensial berikutnya adalah mengestimasi Z (s3 ) berdasarkan data Z (s1 ) dan Z (s 2 ) . Zˆ (s3 ) = ρ13 Z (s1 ) + ρ 23 Z (s 2 )
(2.30)
Estimasi Z berdasarkan tiga data, Z (s1 ) , Z (s 2 ) , dan Z (s3 ) diperoleh dari persamaan :
(
)
(2.31)
− ρ 01 ρ13 − θ 2 ( ρ 23 − ρ12 ρ13 ) )
(2.32)
(3) (2) (2) Zˆ ( s 0 ) = Zˆ ( s 0 ) + θ3 Z ( s3 ) − Zˆ ( s3 )
dengan
θ3 =
1 1 − ρ 212
(ρ
03
Penambahan data baru akan memberikan koreksi terhadap estimasi sebelumnya. Interpolasi terus dilakukan sampai seluruh data yang berada dalam lapangan spasial terikutsertakan dalam estimasi dengan menggunakan persamaan (2.24).
13
II.1.2.4 Kriging untuk Data yang Ditransformasi
Misalkan {Z(s), s ∈ D} adalah suatu proses spasial, dan misalkan {Y(s), s ∈ D} adalah proses spasial intrinsik yang memenuhi Z (s ) = φ ( Y (s ))
(2.33)
dengan φ adalah fungsi yang dapat diturunkan dua kali. Penaksir ordinary kriging bagi Z ( s 0 ) dapat dituliskan sebagai
(
)
ˆ ( s ) + φ '' ( µˆ ) ( σ 2 ( s ) / 2 − m ) Zˆ ( s 0 ) = φ Y 0 Y Y,K 0 Y
(2.34)
ˆ ( s ) adalah penaksir ordinary kriging bagi Y ( s ) , µˆ adalah taksiran dengan Y 0 0 Y nilai mean Y, m Y adalah faktor pengali Lagrange dari sistem persamaan ordinary kriging, dan σ 2Y,K ( s 0 ) adalah variansi kriging di titik s 0 . Sedangkan variansi kriging untuk Z ( s 0 ) dapat dituliskan sebagai 2 σ 2Z,K (s 0 ) = ( φ '(µ Y ) ) σY,K (s 0 ) 2
(2.35)
Uraian tentang hal di atas dapat dilihat pada [4]. Dengan cara yang sama dapat diperoleh penaksir bagi mean Z ( s ) , yaitu ⎛ σ2 mµ ⎞ µˆ Z = φ(µˆ Y ) + φ ''(µˆ Y ) ⎜ Y − Y ⎟ 2 ⎠ ⎝ 2
(2.36)
Dengan σ2Y adalah variansi Y(s), dan mµY faktor pengali Lagrange dari sistem ˆ ( s ) diperoleh melalui sequential persamaan kriging mean. Jika nilai taksiran Y 0 kriging, maka nilai taksiran bagi Z ( s 0 ) dapat dihitung dengan menggunakan persamaan n
(
)
ˆ s ) = φ Y( ˆ s ) + µ − φ ( µ ) − φ '' ( µ ) Z( 0 0 Z Y Y
n
∑∑ α α C(s i =1 j=1
i
j
i
− sj)
2
(2.37)
n −i
dengan α i = θi − ∑ θi + k ρi ,i + k untuk i = 1, 2, …, (n-1) dan α n = θn . k =1
Sedangkan variansi kriging untuk Z ( s 0 ) dapat dituliskan sebagai 2 σ 2Z,K (s 0 ) = ( φ '(µ Y ) ) σY,K (s 0 ) + (−µ Z + φ(µ Y )) 2 2
14
(2.38)
II.2. Kopula Model semivariogram dapat dibentuk dengan menggunakan kopula. Sebelum itu akan dibahas sedikit tentang kopula. Definisikan himpunan bilangan real yang diperluas sebagai R , dengan R = R ∪ {−∞} ∪ {∞} , dan R 2 mendefinisikan ruang dua dimensi bilangan real yang diperluas. Untuk dua vektor x = ( x1, x2 ) dan y = ( y1, y2 ) di R 2 , kita katakan x ≤ y , jika x1 ≤ y1 dan x2 ≤ y2 . Definisi 2.2.1. (Persegi Panjang di R 2 )
Suatu persegi panjang atau interval di R 2 adalah perkalian silang dari dua interval di R 2 , dalam bentuk B = [ x1 ,x2 ] × [ y1 , y2 ] ,
(2.39)
dengan x1 < x2 dan y1 < y2 . Titik-titik sudut persegi panjang B adalah titik-titik ( x1, y1 ), ( x2 , y1 ), ( x1 , y2 ), dan ( x2 , y2 ) . Y
y2 B y1 x1
x2
X
Gambar II.2. Persegi panjang B di R 2
Definisi 2.2.2. (Volume-H)
Misalkan S1 , S2 subhimpunan berupa interval tak kosong di R dan H : R 2 → R adalah suatu fungsi dengan daerah asalnya yaitu Dom( H ) = S1 × S2 . Misalkan
15
B = [ x1 ,x2 ] × [ y1 , y2 ] adalah suatu persegi panjang di Dom( H ) . Maka volume-H
dari persegi panjang B didefinisikan sebagai: VH ( B ) = H ( x2 , y2 ) − H ( x2 , y1 ) − H ( x1 , y2 ) + H ( x1 , y1 )
(2.40)
Jika kita definisikan turunan pertama dari H pada persegi panjang B sebagai ∆ xx2 H ( x, y ) = H ( x2 , y ) − H ( x1 , y ) ,
(2.41)
∆ yy2 H ( x, y ) = H ( x, y2 ) − H ( x, y1 ) .
(2.42)
1
1
Maka volume-H dari persegi panjang B merupakan turunan kedua dari H pada persegi panjang B, yaitu
VH ( B) = ∆ yy2 ∆ xx2 H ( x, y ) . 1
1
(2.43)
Definisi 2.2.3. (2-increasing)
Misalkan H fungsi bernilai real. H dikatakan 2-increasing jika VH ( B) ≥ 0 untuk semua persegi panjang B di R 2 yang semua titik ujungnya berada di Dom(H).
Definisi 2.2.4. (Grounded)
Misalkan S1 , S2 subhimpunan tak kosong dari R dan H : R 2 → R adalah fungsi sedemikian sehingga Dom( H ) = S1 × S2 . Misalkan S1 , S2 mempunyai elemen terkecil masing-masing a dan b . Maka H dikatakan grounded jika
H (a, y ) = 0 = H ( x, b) , untuk setiap x ∈ S1, y ∈ S2 ,
(2.44)
akibatnya jika H grounded, maka
VH ( B) = H ( x, y ) , untuk setiap B = [a, x] × [b, y ] ⊂ Dom( H ) .
(2.45)
karena VH ( B) = H ( x, y ) − H ( x, b) − H (a, y ) + H (a, b) = H ( x, y ) − 0 − 0 + 0 = H ( x, y ) . Definisi 2.2.5. (Subkopula)
Subkopula dua dimensi (2-subcopula) adalah fungsi C ' yang memenuhi beberapa sifat berikut: a. Dom(C ') = S1 × S2 , dengan S1 dan S2 adalah subhimpunan dari I = [0,1] ,
16
b. C ' grounded dan 2-increasing, c. Untuk setiap x ∈ S1 dan y ∈ S2 , berlaku
C '( x,1) = x dan C '(1, y ) = y .
(2.46)
Berdasarkan definisi di atas, 0 ≤ C '(u , v) ≤ 1 untuk setiap (u,v) ∈ Dom(C’ ). Dengan demikian, Range(C ') adalah subhimpunan dari I.
Definisi 2.2.6. (Kopula)
Kopula dua dimensi (2-copula) adalah subkopula 2 dimensi (2-subcopula) C dengan doamin I2.
Definisi di atas setara dengan pernyataan bahwa kopula dua dimensi (2-copula) adalah fungsi C : I 2 → I yang memenuhi sifat-sifat: a. untuk setiap ( x, y ) ∈ I 2 , berlaku
C ( x, 0) = 0 = C (0, y ) ,
(2.47)
C ( x,1) = x dan C (1, y ) = y .
(2.48)
dan
b. untuk setiap u = ( x1 , y1 ), v = ( x2 , y2 ) ∈ I 2 sedemikian sehingga u ≤ v , maka berlaku
C ( x2 , y2 ) − C ( x2 , y1 ) − C ( x1, y2 ) + C ( x1, y1 ) ≥ 0 .
(2.49)
Himpunan dari semua copula dua dimensi didefinisikan sebagai C2. Mengingat, C (u , v) = C (u , v) − C (u , 0) − C (0, v) + C (0, 0) = VC ([ 0, u ] × [ 0, v ])
(2.50)
maka hal ini menunjukkan bahwa C (u , v ) adalah pengaitan suatu bilangan di I terhadap persegi panjang [ 0, u ] × [ 0, v ] .
17
II.2.1 Kopula dan Variabel Acak
Pada subbab ini akan diperkenalkan mengenai teorema Sklar, suatu teorema yang mengaitkan antara copula dengan teori dan aplikasinya dalam ilmu statistika. Namun, sebelumnya akan diperkenalkan beberapa simbol dan definisi yang akan digunakan dalam menjelaskan Teorema Sklar. Peluang bahwa suatu variabel acak Z lebih kecil atau sama dengan z, ditulis P(Z ≤
z) adalah F(z). Nilai F(z) berada di antara 0 dan 1, selanjutnya F(z) disebut dengan fungsi distribusi.
Definisi 2.2.1.1 (Fungsi Distribusi)
Fungsi distribusi (marginal) adalah suatu fungsi F dengan Dom( F ) = R sedemikian sehingga: 1. F fungsi tak turun. 2. F (−∞) = 0 dan F (∞) = 1
Definisi 2.2.1.2 (Fungsi Distribusi Gabungan)
Fungsi distribusi gabungan adalah suatu fungsi H dengan Dom( H ) = R 2 sedemikian sehingga: 1. H fungsi 2-increasing. 2. H ( x, −∞) = H (−∞, y ) = 0 dan H (∞, ∞) = 1 3. H ( x, ∞) = F ( x) dan H (∞, y ) = G ( y ) , dengan F dan G masing-masing adalah fungsi distribusi marginal dari X dan Y. Lema 2.2.1.1 [9]
Misalkan H adalah fungsi distribusi gabungan dengan fungsi distribusi marginalnya masing-masing F dan G, maka terdapat subkopula C ' tunggal sedemikian sehingga: 1. Dom(C ') = Range( F ) × Range(G ) , 2. Untuk setiap x, y ∈ R , H ( x, y ) = C ' ( F ( x), G ( y ) )
18
Bukti:
Karena H memenuhi H ( x, −∞) = 0 = H (−∞, y ) , 2-increasing, dan mempunyai margin F ( x) = H ( x, ∞) dan G ( y ) = H (∞, y ) , maka untuk sebarang titik ( x1 , y1 ) dan ( x2 , y2 ) di R 2 berlaku H ( x2 , y2 ) − H ( x1 , y1 ) ≤ F ( x2 ) − F ( x1 ) + G ( y2 ) − G ( y1 )
Jika F ( x1 ) = F ( x2 ) dan G ( y1 ) = G ( y2 ) , maka H ( x1 , y1 ) = H ( x2 , y2 ) . Selanjutnya, bentuk himpunan pasangan terurut
{( ( F ( x), G( y) ) , H ( x, y) ) | x, y ∈ R} yang mendefinisikan suatu fungsi dua peubah C’ yang bernilai real dengan domain Range( F ) × Range(G ) . Fungsi C’ adalah suatu subkopula berdasarkan sifat-sifat fungsi H. Sebagai contoh, untuk setiap u di Range(F) terdapat suatu x di R
sedemikian
sehingga
F(x)=u.
Maka
C’(u,1)
=
C’(F(x),G( ∞ ))
= H ( x, ∞) = F ( x) = u . Pemeriksaan syarat-syarat lain yang harus dipenuhi suatu subkopula dapat dilakukan dengan cara yang sama. Lema 2.2.1.2 [9]
Misalkan C ' adalah subkopula. Maka terdapat kopula C sedemikian sehingga
C (u , v) = C '(u , v) untuk setiap (u , v) di dom(C ') ,
(2.51)
artinya setiap subkopula dapat diperluas menjadi suatu kopula. Pada umumnya perluasan ini tidak tunggal. Teorema 2.2.1.1 (Teorema Sklar) [9]
Misalkan H adalah fungsi distribusi gabungan dari variable X dan Y, dengan F dan
G masing-masing adalah fungsi distribusi marginal dari X dan Y. Maka terdapat sebuah kopula C sedemikian sehingga untuk setiap x, y ∈ R berlaku H ( x, y ) = C ( F ( x), G ( y ) ) = C (u , v) ,
(2.52)
dengan u = F ( x) dan v = G ( y ) .Jika F dan G kontinu, maka kopula C tunggal. Jika F dan G tidak kontinu, maka kopula C ditentukan secara tunggal pada
Range( F ) × Range(G ) .
19
Sebaliknya, misalkan C adalah sebuah copula, F dan G masing-masing adalah fungsi distribusi dari X dan Y. Maka H pada 2.52 adalah suatu fungsi distribusi gabungan dengan fungsi distribusi marginal F dan G.
Bukti:
Eksistensi kopula pada persamaan 2..52 untuk setiap x, y ∈ R dapat dijelaskan menggunakan lema 2.2.1.1 dan lema 2.2.1.2. Jika F dan G kontinu, berdasarkan lema
2.2.1.1,
maka
terdapat
Dom(C ') = Range( F ) × Range(G ) ,
subkopula karena
F
tunggal
C'
dan
G
dengan
kontinu
maka
Range( F ) = Range(G ) = I atau Dom(C ') = I 2 , ini berarti subkopula tersebut merupakan kopula yang tunggal. Jika F dan G
tidak kontinu, maka terdapat subkopula C ' tunggal
dengan
Dom(C ') = Range( F ) × Range(G ) , maka berdasarkan lema 2.2.1.2 subkopula tersebut dapat diperluas menjadi suatu kopula.
Copula C pada teorema 2.2.1.1 akan dinamakan kopula dari X dan Y, dan dinotasikan
C XY ,
di
mana
kopula
tersebut
dapat
digunakan
untuk
mengidentifikasi dependensi dari variabel acak X dan Y. Salah satu contoh kopula adalah kopula Clayton, yaitu C (u , v) = ( u −α + v −α − 1)
−
1
α
,
dengan α ∈ [ 0, ∞ ) .
II.2.2 Regresi Median
Misalkan X dan Y adalah peubah acak kontinu dengan fungsi distribusi gabungan
H, fungsi distribusi marginal masing-masing F dan G, dan kopula C. Maka U=F(X) dan V=G(Y) adalah peubah acak uniform (0,1) dengan fungsi distribusi gabungan C. Distribusi bersyarat untuk V jika diberikan U = u =F(x) adalah cu (v) = P(V ≤ v | U = u ) = lim
∆u → 0
C (u + ∆u , v) − C (u , v) ∂ = C (u, v) ∆u ∂u
20
(2.53)
Definisi 2.2.2.1 Misalkan X dan Y adalah peubah acak. Untuk x di Ran X , misalkan y = y ( x) adalah solusi persamaan P(Y ≤ y | X = x) =0.5. Maka y = y ( x) adalah persamaan kurva regresi median Y terhadap X.
Perhatikan bahwa
P (Y ≤ y | X = x) = P (V ≤ G ( y ) | U = F ( x)) = P(V ≤ v | U = u ) = cu (v) =
∂ C (u, v) ∂u
Dengan demikian, persamaan kurva regresi median dapat dituliskan sebagai y = G −1 (v ) , dengan v adalah solusi bagi persamaan
∂ C (u, v) = 0.5 . ∂u
II.3 Membangun Model Semivariogram Melalui Regresi Median
Pandang kembali semivariogram eskperimental robust yang diberikan pada persamaan 2.7 yaitu 1 ⎡ ⎧ ⎫⎤ ˆ γ (h) = ⎢ median ⎨ Z (si ) − Z (s j ) 2 : (si , s j ) ∈ N (h) ⎬⎥ ⎩ ⎭⎦ ⎣
4
2 B(h)
1 ⎧ ⎫ Perhatikan bahwa semua anggota ⎨ Z (si ) − Z (s j ) 2 : (si , s j ) ∈ N (h) ⎬ bernilai ⎩ ⎭
positif, dan operasi pemangkatan empat pada bilangan positif tidak mengubah urutan. Dengan demikian , 1 ⎡ ⎧ ⎫⎤ γˆ (h) = ⎢ median ⎨ Z (si ) − Z (s j ) 2 : (si , s j ) ∈ N (h) ⎬⎥ ⎩ ⎭⎦ ⎣ 1 ⎡⎧ ⎫⎤ = median ⎢ ⎨ Z (si ) − Z (s j ) 2 : (si , s j ) ∈ N (h) ⎬⎥ ⎭⎦ ⎣⎩
{
4
4
}
2 B(h)
2 = median ⎡ Z (si ) − Z (s j ) : (si , s j ) ∈ N (h) ⎤ 2 B (h) ⎢⎣ ⎥⎦
⎡ ⎧ Z (s ) − Z (s ) 2 ⎫⎤ ⎪ ⎪ i j ⎢ = median ⎨ : (si , s j ) ∈ N (h) ⎬⎥ ⎢⎪ 2 B(h) ⎪⎭⎥⎦ ⎣⎩
dengan B (h) ≈ 0.457 .
21
2 B(h)
Dari hasil di atas, kita dapat mencari model semivariogram isotropik melalui kurva regresi median variabel
Z (s) − Z (s + h)
2
, dengan h ∈ R 2 dan s, s + h ∈ D,
2 B (h)
terhadap h yang merupakan panjang vektor yang memisahkan lokasi dengan koordinat s dan s + h . Misalkan F adalah fungsi distribusi dari h dan G adalah fungsi
distribusi
dari
Z (s) − Z (s + h) 2 B(h)
2
.
Tuliskan
u = F( h ) ,
dan
⎛ z (s) − z (s + h) 2 ⎞ v = G⎜ ⎟ . Misalkan C(u,v) adalah kopula bagi u dan v. Maka ⎜ ⎟ 2 B ( h ) ⎝ ⎠
model
semivariogram
isotropik
dapat
dirumuskan
sebagai
⎧⎪ G−1(v ), h >0 ∂ γ( h )=⎨ , dengan v adalah solusi bagi persamaan C (u, v) = 0.5 . ∂u ⎪⎩ 0, h =0 Melalui metode ini, tidak perlu lagi memaksakan penggunaan model-model semiovariogram baku, seperti yang selama ini sering dilakukan, sehingga akan mengurangi subjektivitas dalam pemilihan model semivariogram.
22
