MENGHITUNG NILAI PROBABILITAS PADA DISTRIBUSI NORMAL MULTIVARIATE

Menghitung Nilai Probabilitas ... (Jaka Nugraha) MENGHITUNG NILAI PROBABILITAS PADA DISTRIBUSI NORMAL MULTIVARIATE Jaka Nugraha Jurusan Statistika FMIPA Universitas Islam Indonesia e-mail: [email protected] Abstrak Menghitung nilai probabilitas variabel random yang mempunyai distribusi multivariat normal, merupakan salah satu kendala dalam model probit. Persamaan probabilitas dalam model probit tidak berbentuk persamaan tertutup dan harus diselesaikan secara numerik maupun simulasi. Metode Gezn merupakan metode yang paling efektif untuk menghitung nilai probabilitas normal multivariat. Metode ini banyak diimplementasikan dalam beberapa “package” pada program R. Metode ini juga telah diimplementasikan dalam model multinomial probit yang populer disebut dengan simulator Geweke-Hajivassiliou-Keane (GHK). Kata kunci : model multinomial probit, faktor Cholesky, random utiliti, deret Taylor PENDAHULUAN Menghitung nilai probabilitas variabel random yang mempunyai distribusi multivariat normal, merupakan salah satu kendala dalam model probit. Persamaan probabilitas dalam model probit tidak berbentuk persamaan tertutup dan harus diselesaikan secara numerik. Beberapa prosedur non simulasi telah banyak digunakan seperti metode kuadrat untuk mendekati hitungan integral menggunakan fungsi pembobot. Metode ini hanya efektif untuk hitungan integral dimensi rendah. Metode ini dapat digunakan untuk integral dengan dimensi tidak lebih dari empat atau lima ( Geweke, 1996). Metode non simulasi yang lain adalah Clark algoritma, yang telah disampaikan oleh Daganzo dkk (1977). Metode ini sangat tidak akurat sebab akurasi dipengaruhi oleh modelnya (Horowitz dkk, 1982). Untuk menghitung probabilitas dalam model probit dapat juga digunakan metode simulasi. Banyak cara simulasi yang dapat digunakan dalam model probit, secara ringkas dijelaskan oleh Hajivassiliou dkk. (1996) dan menyimpulkan bahwa GewekeHajivassiliou-Keane (GHK) adalah metode yang paling baik paling baik. GHK simulator menghasilkan nilai antara nul s/d satu, dan bersifat unbiased dan konsisten. (Paul C. dkk, 2001)

1

Vol. 3, No. 2, Desember 2007: 114 Permasalahan yang sering muncul dalam statistika adalah menghitung probabilitas normal multivariate. Dalam makalah ini akan dibahas metode menghitung nilai probabilitas variabel random yang mempunyai distribusi multivariat normal menggunakan program R. Pembahasan dimulai dari distribusi normal univariat, distribusi normal bivariate, Distribusi normal

multivariat dan terakhir adalah pembahasan

mengenai pengggunakan Metode GHK dalam model Multinomial Probit.

DISTRIBUSI NORMAL UNIVARIAT Distribusi Normal (Gaussian) yang dinotasikan dengan N(μ ,σ2 ) mempunyai mean μ dan variansi σ2 mempunyai densitas f ( x;  , ) 

  x   2  1 exp  12    dengan μ, σ R 2      

(1)

Selanjutnya, jika variabel random Z berdistribusi normal standar ( mean μ = 0 dan variansi

σ2 = 1). Fungsi distribusi F(z0) yang dinotasikan dengan Φ(z0) adalah 1 2

 ( z0 )  P ( Z  z0 ) 



z0







exp  t 2 / 2 dt

(2)

Untuk mendapatkan nilai (z0) untuk suatu nilai z0, memerlukan tranformasi deret Taylor 

exp(  z 2 / 2)   (1)i i 1

z 2i 2i i!

sehingga 1 2

( z )  FZ ( z )  FZ (0) 



1 1  2 2



 exp  t z

0

 (1)i i 0

2



/ 2 dt

z 2i1 (2i  i)2i i!

(3)

(Paolella, 2006). Invers fungsi (z0), yaitu -1(p), 0 < p < 1, juga dapat dihitung.

DISTRIBUSI NORMAL BIVARIATE Variabel random (X,Y) yang berdistribusi normal bivariate, dengan Var(X) = Var(Y) = 1 dan Korr(X,Y) = , fungsi distribusinya dapat dituliskan dalam bentuk

2

Menghitung Nilai Probabilitas ... (Jaka Nugraha)

 ( x 2  2xy  y 2 )  (b,  )    exp 2 1   2 dydx 2 1   2   1

b1

b2





(4)

Tersa dan Wellan, (1990) menghitung probabilitas normal bivariate dengan mendefinisikan probabilitas dalam bentuk L(h, k ,  ) 

 

1 2 1   2

 h

k

 ( x 2  2 xy  y 2 )  exp  dydx 2 1  2  





(5)

dimana (b,) = L(-b1,-b2, ).

PROBABILITAS NORMAL TRIVARIATE Distribusi normal trivariat dapat dituliskan sebagai

(b, R) 

1

 xT R 1 x  exp     2 dx1dx2dx3 b1

(2 )3 / 2 | R |1 / 2

b2

b3

(6)

dimana b=(b1,b2,b3) dan R = (ij) adalah matrik korelasi. Genz (1992) menyajikan distribusi normal Trivariat berdasarkan distribusi normal bivariat dalam bentuk

(b, R) 

1 (2 )1 / 2

 x2  exp   2  F ( x)dx b1

(7)

dimana  b   x b   x   32  31 21  F ( x)    2 2 211 / 2 , 3 2 311 / 2 , 2 1/ 2 2 1/ 2  ( 1   ) ( 1   ) ( 1   ) ( 1   ) 21 31 21 31   

Terdapat dua metode yang digunakan untuk menghitung F(.) pada persamaan (7), yaitu pertama adalah menggunakan transformasi x=-1(t), sehingga  ( b1 )

(b, R)  

0

F ( 1 (t )dt

(8)

dan selanjutnya itegrasi numerik menggunakan aturan Gauss-Legendre pada interval [0,(b1)]. Yang kedua adalah menggunakan metode Drezner (1992), yang menggunakan modifikasi aturan Gauss-Hermite. Aturan tersebut sangat cocok untuk integral dalam bentuk

1 2





0

x2 exp(  ) f ( x)dx 2

3

Vol. 3, No. 2, Desember 2007: 114 sehingga persamaan (7) ditranformasi kedalam interval [0,]. Misal y=x-b1 dan selanjutnya

(b, R) 



1 (2 )1 / 2

 ( y  b1 ) 2  exp   2  F ( y  b1 )dy 0

 y 2 / 2  yb1  exp( b12 / 2) 0 exp   F ( y  b1 )dy (2 )1 / 2  2  

(9)

jika z=-y maka

(b, R) 





exp( b12 / 2)  exp  z 2 / 2  zb1 F (b1  z )dz (2 )1 / 2 0

(10)

Kedua metode sangat akurat untuk menghitung integral pada daerah terbatas (Genzs,1993). Rumus untuk menghitung nilai probabilitas normal trivariat yang lain didasarkan pada rumus Plackett, yang didasarkan atas derivarif parsial (b, R) exp(  f3 (  21) / 2)  (u3 (  21)) 2  21 2 1   21

dimana

f 3 (r ) 

b3 (1  r 2 )  b1 ( 31  r32 )  b2 ( 32  r31) b12  b22  2rb1b2 dan u ( r )  3 2 2 ((1  r 2 )(1  r 2  31  32  2r3132 )1 / 2 (1  r 2 )

Terdapat dua metode Plackett, yaitu yang pertama menggunakan *  1  21 31    * 2 2 * 21  3132  (1  31 )(1  32 ) dan R*    21 1 32   31 32 1   

Metode Plackett yang ke dua adalah (b, R)  (b, R' ) 

1 2



   1

0



exp(  f3 (  21t ) / 2 21

2 2 1   21 t

(uˆ3 (t ))  31

1 0 b12  b32  2rb1b3 dimana f 2 (r )  dan R'  0 1 2 (1  r ) 0 32

uˆ2 (t ) 

exp(  f 2 ( 31t ) / 2 2 2 1  31 t

 (uˆ2 (t )) dt (11)  

0  32  1 

2 2 b3 (1  31 t )  b1t (  21  3132 )  b3 ( 32  t 2  2131) 2 2 2 2 2 2 2 ((1  31 t )(1  31 t   21 t  32  2t 2 31 2132 )1 / 2

4

Menghitung Nilai Probabilitas ... (Jaka Nugraha)

uˆ3 (t ) 

2 2 b3 (1   21 t )  b1t ( 31   2132 )  b2 ( 32  t 2 31 21) 2 2 2 2 2 2 2 ((1  21 t )(1   21 t  31 t  32  2t 2 31 2132 )1 / 2

Dalam kasus ini, (b,R1) dihitung menggunakan (b1)((b2,b3),32). Metode ke dua ini telah diimplementasikan oleh Drezner (1994).

DISTRIBUSI MULTIVARIATE NORMAL Fungsi densitas multivariate normal dapat disajikan dalam bentuk bm  1  .... exp    t  1 d  m a1 a2 am  2  |  | (2 )

1

F ( a, b ) 

b1 b2

(12)

dimana  = (1, 2,...., m)t dan  adalah matrik covarian mxm yang simetrik positif definet. Pada korelasi konstan, ij =  , untuk ij dan ii = 1, nilai F(.) dapat dihitung dengan akurat menggunakan normal independen. (Tong, 90).

1 F (b)  (2 )1 / 2

 t2  m b  t) ( i )dt  exp  2  1  i 1 

(13)

Secara umum, untuk mengitung integral rangkap, fungsi densitasnya ditranformasi menggunakan faktor cholesky, yaitu  = Cy dimana CCt = . Sehingga t-1 = ytCtC-tC-1Cy = yty dan d = |C|dy = ||1/2dy. Karena a   =Cy  b maka i 1

(ai   cij y j ) j 1

cii

i 1

 yi 

(bi   cij y j ) j 1

cii

untuk i=1,2,..,m

Oleh karena itu F ( a, b) 

1 |  | (2 ) m

b '1



a1

exp( 

bm' ( y1 ,..,ym1 ) y12 b '2 ( y1 ) y2 y2 ) exp(  2 ).... ' exp(  1 )dy am ( y1 ,..,ym1 ) 2 a '2 ( y1 ) 2 2

i 1

i 1

j 1

j 1

(14)

dengan ait ( y1 ,..., yi 1 )  (ai   cij y j ) / cii dan bit ( y1 ,..., yi 1 )  (bi   cij y j ) / cii Untuk masing-masing yi dapat ditransformasi menggunakan yi = -1(zi) , dengan

1 ( y )  (2 )1 / 2

 2  exp  2 d y

5

Vol. 3, No. 2, Desember 2007: 114 yang merupakan fungsi distribusi normal standar univariat. setelah transformasi itu, F(a,b) menjadi F ( a, b)  

e1 e '2 ( z1 )



d1 d2 ( z1 )

...

em ( z1 ,..,zm1 )

dm ( z1 ,..,zm1 )

(15)

dz

dengan i 1 i 1     di ( z1 ,..., zi 1 )  [ai   cij  1 ( z j )] / cii  dan ei ( z1 ,..., zi 1 )  [bi   cij  1 ( z j )] / cii  i 1 i 1    

Bentuk integran ini lebih sederhana dibanding integran aslinya. Integrasi wilayah (region) adalah lebih complek, sehingga tidak bisa ditangani secara langsung menggunakan algoritma numerik standar untuk integral rangkap. Penyelesaian masalah ini menggunakan tranformasi zi = di + wi (ei-di). Setelah tranformasi akhir ini 1

1

1

0

0

0

F (a, b)  (e1  d1 ) (e2  d 2 )...... (em  d m ) dw

(16)

i 1   dengan di  [ai   cij  1 (d j  w j (e j  d j ))] /cii  dan i 1   i 1   ei   (bi   cij  1 (d j  w j (e j  d j ))] /cii  i 1  

Integrasi pertama terhadap variabel w1 , w2 ...dst. 1

1

1

0

0

0

F (a, b)  (e1  d1 ) (em  d m )...... (e2  d 2 )  dw1....dwm

(17)

Integral ini dapat dinyatakan dalam algoritma sebagai berikut 1. Input , , ,a,b dan Nmax 2. menghitung faktor Cholesky bagian bawah, C untuk . 3. nilai awal Itsum=0, N=0, Varsum=0, d1 = (a1/c1,1), e1 = (b1/c1,1) dan f1=e1-d1. 4. Repeat a. Membangkitkan variabel random uniform w1,w2,...,wm-1[0,1] b. untuk i=2,3,....,m i 1   yi-1 = -1(di-1 + wi-1(ei-1-di-1)); di   (ai   ci , j y j ) / cii  i 1   i 1   ei   (bi   ci , j y j ) / cii  dan fi = (ei-di)fi-1 i 1  

6

Menghitung Nilai Probabilitas ... (Jaka Nugraha) c. N=N+1; =(fm – Itsum)/N; Itsum=Itsum+; Varsum=(N-2)Varsum/N + 2 dan Error =  Varsum Until Error <  atau N = Nmax 5. Output : Itsum, Error,dan N Parameter input  adalah faktor confidensi montecarlo untuk error standar. Misalnya untuk  = 2,5, kita mengharapkan error aktual dalah F menjadi lebih kecil dari Error 99%. Nmax adalah parameter input yang merupakan total jumlah perhitungan. Untuk ai = -  atau bi =  maka di = 0 atau ei =1. Dalah beberapa aplikasi, kadang kita menghitung integral

G9 g ) 

1 |  | (2 )

m

b1 b2

bm

a1 a2

am

  ....

 1  exp    t  1 g ( )d  2 

(18)

Untuk menghitung integral ini memerlukan modeifikasi dari algoritma di atas. Karena g() yang tergantung pada m dan wm. Sehingga perlu menghitung  = Cy dan menghitung fm(g()) dalam langkah 4c.

MENGHITUNG NILAI PROBABILITAS DALAM PROGRAM R Terdapat beberapa software yang cukup reliabel dan efisien untuk kasus m=1 dan m=2 [Drezner dkk, 1990] dan [Cox dkk, 1991]. Algoritma Schervish [Schervish , 1984] telah diimplementasikan untuk m kurang dari 8. Algoritma ini menggunakan integrasi secara locally adaptive numerical yang didasarkan pada Simpson’s rule. Algoritma ini membutuhkan waktu yang lama untuk m lebih dari 4 (Genz, 1993). Pendekatan lain adalah dengan metode monte Carlo atau metode subregion adaptive. Metode ini cukup reliabel dan efisien untuk menghitung F(.). [Genz, A., 1992]. Genz (1993) telah melakukan perbandingan beberapa metode untuk menghitung nilai probabilitas normal multivariat. Metode yang dibandingkan adalah acceptance-rejection sampling, metode Deák, metode Genz dan metode Schervish. Disimpulkan bahwa metode Genz adalah yang paling efisien. Pada distribusi normal, untuk menghitung fungsi distribusi maupun fungsi densitas dalam program R, dapat digunakan perintah >pnorm(x, μ, σ) : untuk menghitung nilai fungsi densitas normal

7

Vol. 3, No. 2, Desember 2007: 114 >dnorm(x, μ, σ) : untuk menghitung nilai fungsi distribusi normal Dalam program R juga terdapat beberapa paket yang dapat digunakan untuk menghitung nilai probabilitas multivariat. Perhitungan probabilitas yang didasarkan pada algoritma Genz (1992, 1993)

dapat di akses dalam “fmultivar Package”, “mnormt

Package”, “mvtnorm Package”, “adapt package “.Paket lain yang didasarkan pada metode Schervish's dan metode “Joe” dapat dperoleh pada “mprobit Package”. (R Core Team, 2007)

SIMULASI PROBABILITAS DENGAN METODE GHK DALAM MODEL PROBIT Dalam pemodelan yang didasarkan pada Random Utility, seorang pembuat keputusan dinotasikan dengan i, yang berhadapan dengan pilihan sebanyak J anternatif. Pembuat keputusan mempunyai tingkat utiliti (keuntungan) untuk setiap alternatif. Misalkan Uij untuk j=1,…,J adalah utiliti pembuat keputusan i jika memilih alternatif j. Nilai Uij yang sesungguhya tidak diketahui oleh pengamat (peneliti). Pembuat keputusan memilih alternatif yang mempunyai utiliti terbesar, sehingga memilih alternatif k jika dan hanya jika Uik > Uij j  k. Peneliti tidak mengetahui nilai utiliti untuk pembuat keputusan pada masingmasing alternatif. Peneliti hanya mengamati atribut yang ada untuk masing-masing alternatifnya, dan atribut pembuat keputusan yang dinotasikan dengan x ij. Secara fungsi dapat dinotasikan sebagai Vij= V(xij), j yang biasa dinamakan representative utility. Utiliti dapat dipisahkan menjadi komponen terobservasi (V ij) dan komponen tidak terobservasi (ij), Uij = Vij + ij j,i Probabilitas pembuat keputusan i memilih alternatif k adalah pik = P(Vik + ik > Vij + ij ) jk = =

I(Vik + ik > Vij + ij )(i)di jk I((Vij – Vik) – (ij - ik) < 0)(i)di jk

(19)

dengan I(.) merupakan fungsi indikator dan integral terhadap semua nilai i.

8

Menghitung Nilai Probabilitas ... (Jaka Nugraha) asumsi bahwa vektor  it

Model multinomial probit disusun dari

= (i1,….,

1J)

berdistribusi multivariat normal dengan mean nol dan matrik kovariansi . Densitas untuk i adalah

 ( i ) 

1 (2 )

J /2

||

1/ 2

1 exp[   it  1 i ] 2

(20)

Selanjutnya dapat dilakukan transformasi utiliti ke selisih setiap utiliti terhadap alternatif k, ~ ~ U ijk  U ij  U ik , Vijk  Vij  Vik dan ~ijk   ij   ik .

Sehingga persamaan utilitinya menjadi ~ ~ U ijk  Vijk  ~ijk untuk jk ; ~ikt  (~i1,..., ~iJ ) dimana tanda “ . . .” adalah simbol semua kecuali k, sehingga ~ik adalah matrik (J-1)x1

~ ~ dan ~ik ~N(0,  k ). Matrik kovariansi  k dapat diturunkan dari  dan dapat dicari faktor ~ Choleski, Lk sedemikian hingga Lk Ltk =  k . 0 0 0  c11 0   c21 c22 0 0 0 Lk   c31 c32 c33 0 0   .... ... ... ... ... c  J 1 cJ 2 cJ 3 ... ...

...... ....... ....... ....... c( J 1)1

0  0 0  dimana   0 

Persamaan utiliti dapat disusun sebagai ~ ~ U ik  Vik  Lki Dimana ηti = (η1i, . . . , η(J−1)i) adalah vektor yang berdistribusi normal standard ~ ~ ~ ~ ~ ~ independen, ηji ~ N(0, 1) j ., U ikt  (U i1k ,...,U iJk ) dan Vikt  (Vi1k ,...,ViJk ) Model dapat dituliskan secara eksplisit ~ ~ U i1k  Vi1k  c111i , ~ ~ Ui 2k  Vi 2k  c211i  c212i ~ ~ Ui 3k  Vi 3k  c311i  c322i  c333i

............ ~ ~ UiJk  ViJk  cJ 11i  cJ 22i  cJ 33i  ...  cJJ Ji

9

Vol. 3, No. 2, Desember 2007: 114 Probabilitas individu i memilih alternatif k adalah ~ pik  P(U ijk  0) j  k ~ ~ ~ = P(Ui1k  0,Ui 2k  0,....,UiJk  0 )

~   Vi1k = P1i  c11 

~ ~    ( V  c  )  V i 2 k 21 1 i .P2i  1i  i1k   c c11 22  

 .  

~ ~ ~   (Vi 3k  c311i  c322i )  Vi1k  (Vi 2 k  c211i )   P 3i  1i  dan 2i  ...   c33 c11 c22   ~ ~   (ViJk  cJ 11i  ...  cJ ( J 1)( J 1)i )  Vi1k  P Ji  1i  ,...,  c33 c11  ..., (J -1)i

~  (Vi ( J 1) k  c( J 1)11i  ....  c( J 1)( J  2)( J  2)i )     c( J 1)( J 1) 

(21)

GHK simulator dapat dihitung menggunakan algoritma berikut ini, (Train, 2003) 1. menghitung

~   Vi1k P1  c11 

~    Vi1k      c   11

   

2. mengambil sebuah nilai η1, yang diberi label ηr1 dari distribusi normal terpotong pada ~  V1ki    . Pengambilan dapat dilakukan sebagai berikut :    c  11  (a) mengambil dari distribusi standard uniform μr1. ~   r 1  r  V1ki    (b) menghitung η 1 =  1     c 11    3. menghitung :

~ ~    (Vi 2 k  c211r )   (Vi 2 k  c211 ) r  P2  1  1     c22 c22     Mengambil sebuah nilai η2, diberi label ηr2 dari normal standart terpotong ~ (Vi 2 k  c211r )  . Pengambilan ini dapat dilakukan sebagai berikut : c22 (a) mengambil dari standard uniform μr2.

10

Menghitung Nilai Probabilitas ... (Jaka Nugraha) ~   (V  c211r )    (b) menghitung ηr2 =  1  2r   i 2 k   c 22   

5. menghitung

~ ~    (Vi 3k  c311r  c322r )   (Vi 3k  c311  c322 ) r r  P3  1  1 ,2  2     c33 c33     6. dan seterusnya untuk semua alternatif kecuali k.. 7. probabilitas simulasi untuk pengambilan ke-r dari η1, η2, . . . adalah dihitung sebagai ~ ~ ~   V    (Vi 2 k  c211r )    (Vi 3k  c311r  c322r )    .... pikr   i1k .    c c c 22 33  11      ~   (ViJk  cJ 11r  cJ 22r  ...  cJ ( J 1)(rJ 1) )   ....    c JJ  

8. mengulangi langkah 1-7 beberapa kali, untuk r = 1, . . . , R. 9. probabilitas simulasinya adalah 1 R   pik   pikr R r 1

Algoritma ini sama dengan Algoritma Genz, yaitu dengan batas bawahnya adalah – Inf. ~ ~ pik  P(U ijk  0)  P(  U ijk  0) j  k ~ ~ ~ = P(  Ui1k  0,  Ui 2k  0,....,  UiJk  0 ) ~ ~ ~ = P(  ~i1k  Vi1k ,  ~i 2k  Vi 2k ,....,  ~iJk  ViJk )

sehingga dalam algoritma Genz, diambil nilai d i = 0 dan fi=eifi-1. Sehingga nilai Itsum  sama dengan nilai pik . Misalkan Pengamatan dilakukan dengan mengambil model untuk tiga alternatif dengan memasukkan variabel atribut pembuat keputusan (Xi) dan variabel atribut masing masing alternatif (Zij). Variabel Xi biasa disebut variabel sosio ekonomik/geografi, misalnya penghasilan, jenis kelamin, asal daerah, jumlah anak. Sedangkan variabel Z ij misalkan untuk pilihan penggunaan alat tranportasi (Bus, mobil pribadi, sepeda motor) maka Zij dapat berupa waktu tempuh, biaya. Model utilitinya adalah

11

Vol. 3, No. 2, Desember 2007: 114 Uij = Xi j + Zij + ij untuk i=1,2,...,n dan j=1,2,3. Ui1 =  01 + Xi1 + Zi1 + i1 , Ui2 =  02 + Xi 2 + Zi2+ i2 Ui3 =  03 + Xi3 + Zi3 + i3 Dengan mengambil alternatif pertama sebagai base line, maka model terestimasinya yang merupakan tranformasi selisih Uij-Ui1 adalah U*i21 = ( 02 -  01) + Xi(2- 1) + (Zi2 - Zi1)+ (i2-i1) U*i31 = ( 03 -  01) + Xi(3- 1) + (Zi1 - Zi1) + (i3-i1) U*i11 = 0 Dengan matrik kovariansi (cov(i1,i2, i3)) sama dengan :

  11  12  13        12  22  23     13  23  33  Algoritma GHK tersebut dapat dituliskan dalam program R sebagai berikut : #Program mengitung nilai probailitas dengan metode GHK GHK<-function(b) { b01=b[1]; b02=b[2]; b03=b[3]; b1=b[4]; b2=b[5]; b3=b[6] g=b[7]; s11=b[8];s12=b[9];s13=b[10];s22=b[11] s23=b[12];s33=b[13] V21=(b02-b01) + ( b2-b1)*X+g*(Z[,2]-Z[,1]) V31=(b03-b01) + (b3-b1)*X+g*(Z[,3]-Z[,1]) V32=V31-V21 S=matrix(c(s11,s12,s13,s12,s22,s23,s13,s23,s33),3,3) M1= matrix(c(-1,-1,1,0,0,1),2,3) M2=matrix(c(1,0,-1,-1,0,1),2,3) S1= M2%*%S%*%t(M2) S2=M2%*%S%*%t(M2) F1=matrix(,N,1); F2= matrix(,N,1);F3= matrix(,N,1) Prob=matrix(,N,3) for (i in 1:N) { F1= sadmvn(-Inf, c(-V21[i],-V31[i]), mean=rep(0,2),S1) F2= sadmvn(-Inf, c(V21[i],-V32[i]), mean=rep(0,2),S2) F3= 1 – F1-F2 Prob[i,]=c(F1,F2,F3) } Prob }

Misalkan diambil data dengan N= 500, Xi ~ NID(0,1) , Zij ~ NID(0,1) dan i ~ N(0,)

 1 0.2 0.7    =  0.2 1 0.7   0.7 0.7 1   

12

Menghitung Nilai Probabilitas ... (Jaka Nugraha) parameter 01 =2 , 02 =1, 03 =0.5, 1=-2, 2=-1,  3 =-1 dan

=0.8, maka nilai

probabilitas masing-masing alternatif(pilihan) pada masing-masing individu dapat dihitung sebagai berikut: > library(mnormt) >X=rnorm(N,3); prop.pilihan= matrix(,N,3); Z=matrix(,N,3); V=matrix(,N,3) >for (i in 1:3){Z[,i]=rnorm(N)} >prop.pilihan=GHK(c(2,1,0.5,-2,-1,-1,0.8,1,0.2,0.7,1,0.7,1))

KESIMPULAN Seiring dengan kemajuan teknologi komputasi, saat ini memungkinkan menghitung nilai probabilitas dari variabel random yang berdistribusi normal multivariate. Metode yang banyak diaplikasikan dalam program R adalah Metode Genz yang disusun berdasarkan simulasi Monte Carlo. Aplikasi metode ini dalam bidang ekonometrika, adalah untuk simulasi nilai probabilitas dikenal dengan nama metode GHK. Dengan adanya simulasi probabilitas (GHK) ini, maka dengan program R adalah sangat memungkinkan untuk menyusun model multinomial probit.

DAFTAR PUSTAKA Contoyannis P, 2001, An Introduction to Simulation-Based Estimation, Working Paper Department of Economics and Related Studies, University of York. Cox, D.R. dan Wermuth, N. (1991), A Simple Approximationfor Bivariateand Normal Integrals, International Statistics Review 59, pp. 263-269 Daganzo, C., F. Bouthelier, and Y. Sheffi (1977), „Multinomial probit and qualitative choice: A computationally efficient algorithm‟, Transportation Science 11, 338– 358. Drezner, Z. (1992), Computation of Multivariate Normal Integral, ACM Transactions on Mathematics Software 18, pp. 450-460 Drezner, Z. (1994), Computation of the Trivariate Normal Integral, Mathematics of Computation 62, pp. 289-294

13

Vol. 3, No. 2, Desember 2007: 114 Drezner, Z. dan Wesolowsky G.O., (1989), On the Computation of Bivariate Normal Integral, Journal of Statist. Comput. Simul. 35. pp. 101-107 Genz, A. (1992). Numerical Computation of Multivariate Normal Probabilities. J. Computational and Graphical Statist., 1, 141-149. Genz, A. (1993). Comparison of methods for the computation of multivariate normal probabilities. Computing Science and Statistics, 25, 400-405. Geweke, J. (1996), „Monte Carlo simulation and numerical integration‟, in D. Kendrick and J. Rust, eds., Handbook of Computational Economics, Elsevier Science, Amsterdam, pp. 731–800. Hajivassiliou, V., D. McFadden, and P. Ruud (1996), „Simulation of multivariate normal rectangle probabilities and their derivatives: Theoretical and computational results‟,Journal of Econometrics 72, 85–134. Horowitz, J., J. Sparmann, and C. Daganzo (1982), „An investigation of the accuracy of the Clark approximation for the multinomial probit model‟, Transportation Science 16, 382–401 Paolella, Marc S. (2006), Fundamental Probability (A Computational Approach), John Wiley & Sons Ltd, Schervish, M.J. (1984). Multivariate normal probabilities with error bound. Appl. Statist., 33, 81-94. Tong, Y.L. (190), The Multivariate Normal Distribution, Springer-Verlag, New York Train, Kenneth (2003), Discrete Choice Methods with Simulation, UK Press, Cambridge.

14