RANCANGAN D-OPTIMAL MODEL MICHAELIS MENTEN DAN EMAX DENGAN MATLAB
Tatik Widiharih 1 , Sri Haryatmi 2 , Gunardi 3 , Yuciana Wilandari 4
1 Jurusan Statistika FSM UNDIP, widiharih@gmail.com
2 Jurusan Matematika FMIPA UGM, s_kartiko@yahoo.com
2 Jurusan Matematika FMIPA UGM, gunardi@ugm.ac.id
1 Jurusan Statistika FSM UNDIP, yuciana.wilandari@gmail.com

## Abstract
Michaelis Menten and Emax  models    are    widely used  in chemistry,  pharmacokinetics  and pharmacodynamics areas. D-optimal criteria is criteria with the purpuse to minimize the variance of the estimator  of  parameters in  the model.  In  this paper  will  discuss  the D-optimal  design  for Michaelis Menten and Emax models  with    homoscedastics    error  assumtion.    Determination  of D-optimal designs based on Generalied Equivalence Theorem Kiefer-Wolvowitz.  We used minimally supported design with the proportion of  each design point is uniform, lower bound of design region is design point and the others are interior points.
Keywords : D-optimal, Michaelis Menten, Emax, Minimally Supported Design, Homoscedastics

## 1. Pendahuluan
Sebelum melakukan penelitian, peneliti membuat rancangan yang sesuai dengan permasalahan yang dihadapi. Masalah yang sering muncul adalah bagaimana menentukan titik-titik rancangan yang akan dicobakan sehingga memenuhi kriteria yang telah ditetapkan oleh peneliti. Salah satu kriteria tersebut adalah kriteria optimal. Optimal bisa berarti mempunyai galat minimum, biaya minimum, bias minimum, variansi minimum dan sebagainya sesuai dengan permasalahan dan tuuan yang telah ditetapkan.
Kriteria D-optimal merupakan kriteria dengan tujuan meminimalkan variansi dari penduga parameter dalam model. Bila  variansi dari penduga parameter  kecil diharapkan model akan signifikan (hipotesis yang menyatakan nilai parameter sama dengan nol ditolak).  Pada awalnya kriteria D-optimal digunakan untuk model regresi linear. Selanjutnya dikembangkan untuk model nonlinear dengan menggunakan ekspansi deret Taylor [1],[9] . Model nonlinear banyak diaplikasikan pada berbagai bidang terutama untuk memodelkan kurva pertumbuhan.
Model regresi nonlinear yang banyak digunakan diantaranya model sigmoid. Rancangan D-optimal untuk model nonlinear telah banyak dilakukan diantaranya Dette and Pepelyshev [2] menggunakan model sigmoidal untuk empat macam model yaitu  model regresi eksponensial dengan bentuk ε β λ + − = − x e y α ,  model regresi  Weibull dengan bentuk ε β γ λ x e y − − = α , model regresi logistic dengan bentuk ε β α λ + + = − x e y α , dan model regresi Richards dengan bentuk ( ) ε β α γ λ + + = − x e y α .  Li [5] menggunakan model regresi
Gompertz dengan bentuk ( ) ) ( exp . τ γ β − − − = x e y , yang diaplikasikan pada bidang biologi. Model  Michaelis Menten dan Emax banyak diaplikasikan pada bidang kimia, pharmakokinetik dan pharmakodinamik.  Jericevic and Kuster [8] menggunakan model Michaelis Menten untuk menggambarkan hubungan antara konsentrasi dari  substrate dan
kecepatan reaksi. Holford [7] menggunakan model EMAX untuk menggambarkan hubungan antara konsentrasi dari   receptor dan pengaruh dari sejenis obat. Kurva dari model Michaelis Menten dan Emax berbentuk huruf S.
Rancangan D-optimal untuk model kurva pertumbuhan berbentuk huruf S juga telah diteliti diantaranya Li and Majumdar [4] menggunakan fungsi Logistik dengan tiga dan empat parameter, dan rancangan yang digunakan merupakan rancangan minimal. Fang and Hedayat [3] menggunakan model Emax dengan satu komponen melalui pendekatan analisis fungsional. Selanjutnya Li [5] menggunakan fungsi Gompertz  dengan dua dan tiga parameter. Rancangan yang digunakan adalah rancangan minimal dan semua titik rancangan  berupa titik interior dari daerah rancangan. Rancangan D-optimal untuk model kurva pertumbuhan tumor telah diteliti oleh Li and Balakrishnan [6] . Model yang digunakan dikenal dengan model LINEX  dan REGROWTH. Rancangan yang digunakan adalah rancangan minimal dengan titik-titik rancangan semua merupakan titik interior dari daerah rancangan.
Terinspirasi untuk kurva pertumbuhan berbentuk huruf S, sebagai alternatif digunakan fungsi Morgan Mercer Flodin (MMF). Pada model ini penulis telah melakukan penelitian untuk model MMF tanpa intercept dengan pendekatan konsep sistem Tchebysheff [10] . Selanjutnya penelitian  rancangan D-optimal untuk model MMF tiga parameter telah juga dilakukan oleh Wdiharih, et al. [11] .
Pada makalah ini model Michaelis Menten yang akan  digunakan mempunyai bentuk sebagai berikut:
[ ] b x x x y , 0 , 0 , , 3 2 2 3 ∈ > + + = θ θ ε θ θ (1)
Model Emax  yang digunakan mempunyai bentuk sebagai berikut:
[ ] b x x x y , 0 , 1 , 0 , , 4 3 2 2 3 4 4 ∈ ≥ > + + = θ θ θ ε θ θ θ θ (2)
dengan  asumsi    merupakan variabel acak  dengan rata-rata nol dan variansi konstan . Berdasarkan model pada Persamaan (1) dan (2) akan ditentukan titik-titik rancangan D-optimal dengan program Matlab yang berupa M-file untuk mendefinisikan fungsi sasaran, menentukan titik-titik rancangan, menentukan matriks informasi dan inversnya serta untuk mengkonstruksikan grafik dari fungsi variansi terstandar.

## 2. Deskripsi Teoritis.
Model  nonlinear yang digunakan mempunyai bentuk ( ) ε β η + = , x y . Rancangan dengan  p  buah titik rancangan dinotasikan:
      = p p w w w x x x ... ... 2 1 2 1 ξ (3)
dengan: n r w i i = i r : banyaknya ulangan untuk titik rancangan x i
n : banyaknya pengamatan  dan
∑ = = p i i r n 1
1 1 = ∑ = p i i w
Matriks informasi dari rancangan (3) adalah :
∑ = = d i T i x h x h w M 1 ) ( ) ( ) ( ξ (4)
merupakan matriks simetri berukuran kxk dengan k adalah banyaknya parameter dalam model,  dengan . ) , ( ,..., ) , ( , ) , ( ) ( 2 1 T k x x x x h       ∂ ∂ ∂ ∂ ∂ ∂ = β β η β β η β β η Matriks  informasi pada
Persamaan (4) masih mengandung parameter. Dalam melakukan percobaan, peneliti telah mempunyai informasi tentang nilai parameter, sehingga bila nilai parameter dimasukkan pada Persamaan (4), determinan dari matriks informasi merupakan fungsi dengan peubah ganda dalam x ( x 1 , x 2 ,...,x p ). Fungsi dispersi (variansi terstandar) yang bersesuaian dengan rancangan ξ adalah:
) ( ) ( ) ( ) , ( 1 x h M x h x d T ξ ξ − = (5)
Nilai  dari  variansi terstandar  untuk  rancangan D-optimal antara nol sampai dengan k ( k adalah banyaknya parameter dalam model) dan maksimum dari variansi terstandar terjadi pada titik-titik rancangan. Berikut definisi dari rancangan dengan kriteria D-optimal.
Definisi 1. Atkinson et al [1]
Kriteria D-optimal adalah kriteria yang diperoleh dengan cara meminimumkan { } ) ( ξ Ψ M dengan: { } ) ( log ) ( log ) ( 1 ξ − = ξ = ξ Ψ − M M M
Berdasarkan Definisi 1, kriteria D - optimal diperoleh dengan cara memaksimalkan ) ( ξ M . Untuk menunjukan bahwa rancangan ξ seperti pada Persamaan (3) memenuhi kriteria D-optimal digunakan  Teorema  Equivalensi  Tergeneral dari Kiefer and Wolfowitz [9] berikut.
Teorema 1. Teorema Equivalensi Tergeneral Kiefer and Wolfowitz [9]
Tiga kondisi berikut adalah equivalen:
1. Rancangan * ξ memaksimalkan ) ( ξ M
2. Rancangan * ξ meminimumkan ( ) ξ . max x d x 3. ( ) k x d x = ξ . max , dengan k adalah banyaknya parameter dalam model.
Secara ringkas Teorema Equivalensi Tergeneral dapat dinyakan sebagi berikut:
Rancangan * ξ memaksimalkan k x h M x h x d M T ≤ ξ = ξ ⇔ ξ − ) ( ) ( ) ( ) , ( ( 1 (6)
Untuk menunjukkan bahwa rancangan    seperti pada Persamaan (3) merupakan rancangan D-optimal dengan melakukan cheking bahwa terpenuhinya:
k x h M x h x d T ≤ ξ = ξ − ) ( ) ( ) ( ) , ( 1
3. Hasil dan Pembahasan 3.1 Rancangan D-optimal Model Michaelis Menten.
Model Michaelis Menten yang digunakan adalah:
[ ] b x x x y , 0 , 0 , , 3 2 2 3 ∈ > + + = θ θ ε θ θ x x x + = 2 3 ) , ( θ θ θ η Turunan parsial dari ) , ( θ η x terhadap parameternya adalah:
( ) ( ) T x x x x x h       + + − = 2 2 2 3 , ) ( θ θ θ
Grafik dari model Michaelis Menten  untuk 1 3 = θ pada  beberapa nilai 2 θ seperti pada Gambar 1.
Gambar 1 . Model Michaelis Menten   untuk 1 3 = θ pada  Beberapa Nilai 2 θ
Untuk menentukan rancangan D-optimal (menentukan titik-titik rancangan dan proporsinya sehingga terpenuhi kriteria D-optimal ) dari model seperti Persamaan (1) dengan menggunakan Teorema berikut:
Teorema 2. (Widiharih, et al . [11] ).
Rancangan D-optimal model Michaelis Menten  merupakan rancangan minimal dengan titik-titik rancangan adalah: b x b b x = + = 2 2 2 1 , 2 θ θ dengan proporsi yang sama
yaitu 2 1 .
Berdasarkan Teorema 2, karena menggunakan rancangan minimal maka:
        = 2 1 2 1 1 b x ξ
Elemen matriks informasi adalah:
( ) ∑ = + θ θ = 2 1 4 2 2 2 3 11 2 1 i i i x x m ( ) ∑ = + θ = 2 1 2 2 2 22 2 1 i i i x x m ( ) ∑ = + θ θ − = 2 1 3 2 2 3 12 2 1 i i i x x m b x = 2
Berdasarkan Teorema 1, dapat dibuat program Matlab untuk menenitukan titik-titk rancangan D-optimal dan  cheking Teorema Equivalensi Tergeneral (Lampiran 1). Pada kasus 3976 , 25 , 3262 , 2 3 2 = = θ θ dengan daerah rancangan [0 , 1.5], diperoleh  titik-titik rancangan 5671 , 0 1 = x dan 5 , 1 2 = x .
Matriks informasi, invers matriks informasi dan fungsi dispersi  (variansi terstandar) adalah:
      − − = 0961 , 0 6787 , 0 6787 , 0 8659 , 4 ) ( ξ M       = − 6986 , 719 3838 , 100 3838 , 100 2071 , 14 ) ( 1 ξ M ) ( ) ( ) ( ) , ( 1 x h M x h x d T ξ ξ − =             + + −             + + − = ) 3262 , 2 ( ) 3262 , 2 ( 3976 , 25
6986 ,
719 3838 , 100 3838 , 100 2071 , 14 ) 3262 , 2 ( ) 3262 , 2 ( 3976 , 25 2 2 x x x x x x x x
Grafik dari fungsi variansi terstandar seperti pada Gambar 2 berikut:
Gambar 2 . Grafik ) , ( ξ x d untuk Model (1) dengan
3976 , 25 , 3262 , 2 3 2 = = θ θ pada Daerah Rancangan [0, 1.5]
Berdasar Gambar 2, terlihat bahwa
2 ) , ( ≤ ξ x d sehingga memenuhi Teorema Equivalensi Tergeneral, dapat disimpulkan titik-titik rancangan 5671 , 0 1 = x dan 5 , 1 2 = x merupakan titik-titik rancangan D-optimal model (1)  dengan 3976 , 25 , 3262 , 2 3 2 = = θ θ pada daerah rancangan [0, 1.5]  .

## 3.2 Rancangan D-optimal Model Emax.
Model Emax yang digunakan adalah:
[ ] b x x x y , 0 , 1 , 0 , , 4 3 2 2 3 4 4 ∈ ≥ θ > θ θ ε + + θ θ = θ θ 4 4 2 3 ) , ( θ θ + θ θ = θ η x x x (7)
Turunan parsial dari ) , ( θ η x terhadap parameternya adalah: ( ) ( ) ( ) T x x x x x x x x h         + + + − = 2 2 3 2 2 2 2 3 4 4 4 4 4 4 ) ln( , , ) , ( θ θ θ θ θ θ θ θ θ θ θ θ θ Grafik dari model Emax  untuk 3262 , 2 , 1385 , 10 4 2 = = θ θ pada  beberapa nilai 3 θ seperti pada Gambar 3.
Gambar 3 . Model Emax  untuk
3262 , 2 , 1385 , 10 4 2 = = θ θ pada  Beberapa Nilai 3 θ
Untuk  menentukan titik-titik rancangan D-optimal model EMAX menggunakan teorema berikut:
Teorema 3. (Widiharih, et al , [11] ).
Rancangan D-optimal model EMAX merupakan rancangan minimal dengan titik-titik rancangan adalah x 1 , x 2 dan x 3 = b dengan proporsi sama yaitu 3 1 , x 1 , x 2 ,  merupakan nilai yang memaksimumkan ( ) ( ) ( ) B A x x x M + ∝ 4 2 3 2 1 , θ θ ξ , dengan: ( ) ( ) ( ) ( ) ( ) ( ) ∑ = + + + − = 3 1 2 2 4 2 4 2 2 4 4 4 ln ln i k j i j i x x x x x A θ θ θ θ θ θ ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ∑ = + + + − − + = 3 1 3 2 3 2 4 2 2 4 4 4 ln ln ln ln ln ln ln i k j i k j i k i j i x x x x x x x x x x B θ θ θ θ θ θ i ≠ j ≠ k , j , k = 1, 2, 3, x 3 = b
Berdasar Teorema  2 dapat diperoleh:
        = 3 1 3 1 3 1 2 1 b x x ξ Elemen matriks informasi adalah: ( ) ∑ = + = 3 1 4 2 2 3 11 4 4 3 1 i i i x x m θ θ θ θ ( ) ∑ = + = 3 1 2 2 22 4 4 3 1 i i i x x m θ θ θ ( ) ( ) ∑ = + = 3 1 4 2 2 2 2 3 2 2 33 4 4 ln 3 1 i i i i x x x m θ θ θ θ θ ( ) ∑ = + − = 3 1 3 2 2 3 12 4 4 3 1 i i i x x m θ θ θ θ ( ) ( ) ∑ = + = 3 1 4 2 2 2 3 2 13 4 4 ln 3 1 i i i i x x x m θ θ θ θ θ ( ) ( ) ∑ = + = 3 1 3 2 2 2 3 2 23 4 4 ln 3 1 i i i i x x x m θ θ θ θ θ
dengan x 3 = b
Berdasarkan Teorema 3, dapat dibuat program Matlab untuk menentukan titik-titk rancangan D-optimal dan cheking Teorema Equivalensi Tergeneral (Lampiran 2). Pada kasus 008629 , 1 , 77242 , 99 , 1385 , 10 4 3 2 = = = θ θ θ , diperoleh diperoleh titik-titik  rancangan 1451 , 2 1 = x , 9973 , 12 2 = x dan 0000 , 50 3 = x
Matriks informasi , invers matriks informasi dan fungsi dispersi (variansi terstandar) adalah:
          − − − − = 3095 , 2 0274 , 0 0799 , 0 0274 , 0 0004 , 0 0009 , 0 0799 , 0 0009 , 0 0032 , 0 10 ) ( 3 ξ M           − − − − = − 0121 , 0 6350 , 0 1184 , 0 6350 , 0 4450 , 44 0995 , 3 1184 , 0 0995 , 3 3606 , 2 ) ( 1 ξ M ) ( ) ( ) ( ) , ( 1 x h M x h x d T ξ ξ − =         
         + + + −           − − − −       + + + − = 2 2 3 2 2 2 2 3 2 2 3 2 2 2 2 3 ) ( ) ln( ) ( ) ( 0121 , 0 6350 , 0 1184 , 0 6350 , 0 4450 , 44 0995 , 3 1184 , 0 0995 , 3 3606 , 2 ) ( ) ln( ) ( ) ( 4 4 4 4 4 4 4 4 4 4 4 4 θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ x x x x x x x x x x x x x x dengan: 008629 , 1 , 77242 , 99 , 1385 , 10 4 3 2 = = = θ θ θ
Grafik dari fungsi variansi terstandar seperti pada Gambar 4 berikut:
Gambar 4 . Grafik ) , ( ξ x d untuk Model (2) dengan 008629 , 1 , 77242 , 99 , 1385 , 10 4 3 2 = = = θ θ θ
pada Daerah Rancangan [0, 50]

## 4. Kesimpulan
Berdasarkan Teorema 2  dan  3, penentuan rancangan D-optimal model Michaelis Menten dan Emax mempunyai bentuk yang rumit. Determinan matriks  informasi mempunyai bentuk yang kompleks, penentuan titik maksimum dari determinan matriks informasi tidak bisa dilakukan secara langsung. Pendekatan numeris sebagai alternatif dari masalah optimasi, dalam hal ini digunakan program Matlab untuk memudahkan  para pengguna.

## DAFTAR PUSTAKA
1. Atkinson, A.C., Donev, A.N. and Tobias, R.D., Optimum Experimental Designs, with SAS, OXFORD University Press, 2007.
2. Dette, H and Pepelyshev, A., Efficient Experimental Designs for Sigmoidal Growth Models. Journal of Statistical Planning and  Inference , 2008, 138: 2-17.
3. Fang, X., and Hedayat, A.S., Locally D-optimal Designs Based on a Class of Compossed Models Resulted from Blending Emax and One-compartment Models, The Annals of Statistics , 2008, Vol.36, No.1: 428-444.
4. Li,G., and Majumdar, D., D-optimal Designs for Logistic Models with Three and Four Parameters, Journal of Statistical Planning and Inference, 2008, 138: 190-1959, .
5. Li, G., Optimal and Eficient Designs for Gompertz Regression Models. Ann Inst Stat Math, DOI 10.1007/s10463-011-03040-y, 2011
6. Li, G., and Balakrishnan, N.,  Optimal Designs for Tumor Regrowth Models, Journal of Statistical Planning and Inference , 2011, 141: 644-654.
7. Holford, N., Pharmacodynamic Principles and the Time Course of Immediate Drug Effect, Department Pharmacology and Clinical Pharmacology, University of Auckland, New Zeland, 2013 .
8. Jericevic, Z.,and Kuster, Z., Nonlinear Optimization of Parameters in Michaelis Menten Kinetics, Croatia Chemica Acta CCACCAA , 2005, Vol. 78, No. 4: 519-523. 9. Kiefer, J., and Wolfowitz, J., The Equivalence of Two Extremum Problems, Can. Jnl. Math, 1960, Vol. 12: 363-36.
10. Scwhabe, R., Optimal Design for Linear and Nonlinear Models, A Short Course. Otto von Guerike University Magdeburg, July 2008, 21-25 .
11. Widiharih,T., Haryatmi,S., dan Gunardi., D-optimal designs for Morgan Mercer Flodin (MMF) Models without Intercept , International Journal of Applied Mathematics and Statistics , 2015a, Vol.53, No. 5:163-171.
12. Widiharih,T., Haryatmi,S., dan Gunardi., D-optimal designs for Morgan Mercer Flodin (MMF) Models With Three Parameters, Proceeding AIP, SEAMS UGM 2015b.
Lampiran 1 . Program Matlab Untuk Menentukan Rancangan D-optimal Model Michaelis Menten.
Model Michaelis Menten yang digunakan dalam program :

## M-FILE UNTUK MENENTUKAN TITIK-TITIK RANCANGAN,MATRIKS INFORMASI DAN
INVERS MATRIKS INFORMASI function mm syms a b lb  x y xx n11 n22 n12 a= ;  % masukkan nilai a b= ; % masukkan nilai b lb= ; % masukkan nilai batas atas daerah rancangan x=(a*lb)/(2*a+lb); ds1=x; y=lb; ds2=y; n11=0.5*[((b*x)^2)/((a+x)^4)+ ((b*lb)^2)/((a+lb)^4)]; m11=simplify(n11); n22=0.5*[x^2/((a+x)^2)+ lb^2/((a+lb)^2)]; m22=simplify(n22); n12=-0.5*[(b*x^2)/((a+x(1))^3)+ (b*lb^2)/((a+lb)^3)]; m12=simplify(n12); M=[m11 m12;m12 m22]; inves=inv(M) save in inves
M-FILE UNTUK MENGGAMBAR VARIANSI TERSTANDAR function gmm syms x invs=load('in.mat'); d1=struct2cell(invs); d2=cell2mat(d1); a= ; % masukkan nilai a b= ; % masukkan nilai b a1=-b*x/(a+x)^2; a11=simplify(a1); a2=x/(a+x); a22=simplify(a2); d=[a1;a2]'*d2*[a1;a2]; ezplot(d) axis([bbx bax bby bay]) grid on
% bbx : batas bawah dari x % bax : batas atas dari x % bby : batas bawah dari y % bay : batas atas dari y
Lampiran 2 . Program Matlab Untuk Menentukan Rancangan D-optimal Model Emax.
Model Emax yang digunakan dalam program :
M-FILE UNTUK MENDEFINISIKAN DETERMINAN MATRIKS INFORMASI function f=emax(x,b,c,d) b= ; %masukkan nilai b c= ; %masukkan nilai c d= ; %masukkan nilai d x1=x(1); x2=x(2); x3=x(3); m11=0.333*c^2*[x1^(2*d)/(b+x1^d)^4+x2^(2*d)/(b+x2^d)^4 +x3^(2*d)/(b+x3^d)^4]; m22=0.333*[x1^(2*d)/(b+x1^d)^2+x2^(2*d)/(b+x2^d)^2 +x3^(2*d)/(b+x3^d)^2]; m33=0.333*b^2*c^2*[x1^(2*d)*[log(x1)]^2/(b+x1^d)^4
+x2^(2*d)*[log(x2)]^2/(b+x2^d)^4
+x3^(2*d)*[log(x3)]^2/(b+x3^d)^4]; m12=-0.333*c*[x1^(2*d)/(b+x1^d)^3 +x2^(2*d)/(b+x2^d)^3+x3^(2*d)/(b+x3^d)^3]; m13=-0.333*b*c^2*[x1^(2*d)*[log(x1)]/(b+x1^d)^4 +x2^(2*d)*[log(x2)]/(b+x2^d)^4 +x3^(2*d)*[log(x3)]/(b+x3^d)^4]; m23=0.333*b*c*[x1^(2*d)*[log(x1)]/(b+x1^d)^3 +x2^(2*d)*[log(x2)]/(b+x2^d)^3 +x3^(2*d)*[log(x3)]/(b+x3^d)^3]; M=[m11 m12 m13; m12 m22 m23; m13 m23 m33]; f=-det(M)

## M-FILE UNTUK MENENTUKAN TITIK-TITIK RANCANGAN,MATRIKS INFORMASI DAN
INVERS MATRIKS INFORMASI function x=maksemax(Lb,Ub,b,c,d) Lb=[ ; ; ];  %masukkan batas bawah untuk x1,x2 dan x3 Ub=[ ; ; ];  % masukkan batas atas untuk x1,x2 dan x3 x0=[ ; ; ];  % masukkan nilai awal untuk x1,x2 dan x3 [x,nilaimaximum,ExitFlag,output]=fmincon('emax',x0,[],[],
[],[],Lb,Ub)
b= ;
%masukkan nilai b c= ; %masukkan nilai c
d= ; %masukkan nilai d x1=x(1); x2=x(2); x3=x(3); m11=0.33333*c^2*[x1^(2*d)/(b+x1^d)^4+x2^(2*d)/(b+x2^d)^4 +x3^(2*d)/(b+x3^d)^4]; m22=0.33333*[x1^(2*d)/(b+x1^d)^2+x2^(2*d)/(b+x2^d)^2 +x3^(2*d)/(b+x3^d)^2]; m33=0.33333*b^2*c^2*[x1^(2*d)*[log(x1)]^2/(b+x1^d)^4 +x2^(2*d)*[log(x2)]^2/(b+x2^d)^4 +x3^(2*d)*[log(x3)]^2/(b+x3^d)^4]; m12=-0.33333*c*[x1^(2*d)/(b+x1^d)^3 +x2^(2*d)/(b+x2^d)^3+x3^(2*d)/(b+x3^d)^3]; m13=-0.33333*b*c^2*[x1^(2*d)*[log(x1)]/(b+x1^d)^4 +x2^(2*d)*[log(x2)]/(b+x2^d)^4 +x3^(2*d)*[log(x3)]/(b+x3^d)^4]; m23=0.33333*b*c*[x1^(2*d)*[log(x1)]/(b+x1^d)^3 +x2^(2*d)*[log(x2)]/(b+x2^d)^3 +x3^(2*d)*[log(x3)]/(b+x3^d)^3]; M=[m11 m12 m13; m12 m22 m23; m13 m23 m33]
inves=inv(M) save in inves
M-FILE UNTUK MENGGAMBAR VARIANSI TERSTANDAR function gemax syms x inves=load('in.mat'); d1=struct2cell(inves); d2=cell2mat(d1); b= ; %masukkan nilai b c= ; %masukkan nilai c d= ; %masukkan nilai d a1=-c*(x^d)/(b+x^d)^2; a2=(x^d)/(b+x^d); a3=(b*c*x^d*log(x))/(b+x^d)^2; dd=[a1 a2 a3]*d2*[a1;a2;a3]; ezplot(dd, [bbx,bax]) grid on
% bbx : batas bawah dari x % bax : batas atas dari x