Newton Rhapson Method

Metode Newton Rhapson adalah sebuah metode numerik yang digunakan untuk membantu estimasi parameter sebuah distribusi. Akan dicobakan kedalam 4 distribusi yaitu Distribusi Gamma, Distribusi Weibull,Distribusi Normal dan distribusi eksponensial.

1. Distribusi Gamma

Fungsi pdf Gamma adalah \(f(x)\)=\(\frac{1}{\Gamma (\alpha) \beta ^{\alpha}} x^{\alpha - 1} e^{- \frac{x}{\beta}}\)

nr.gamma <- function(x,theta0){
  n=length(x)
  #converg<-1000
  e=10^-6
  alpha_est=c()
  beta_est=c()
  error=c()
  #estimasi parameter gamma
  alpha0 <- theta0[1]
  beta0 <- theta0[2]
  alpha_est[1]<-alpha0
  beta_est[1]<- beta0
  i=1
  # Newton-Raphson hanya diturunkan terhadap alpha
  repeat{
    alpha <- alpha_est[i]
    der1<- -n*log(mean(x)/alpha)-n*digamma(alpha)+sum(log(x))
    der2<- n/alpha-n*trigamma(alpha)
    alpha_est[i+1]<-alpha_est[i]-der1/der2
    error[i+1]<-abs(alpha_est[i+1]-alpha_est[i])
    if (error[i+1]<e){break}
    beta_est <-c(beta0,mean(x)/alpha_est)
    i=i+1
  }
  par(mfrow=c(2,1))
  return(list(cbind(alpha_est,beta_est,error),alpha=alpha,beta=beta_est[length(beta_est)],plot(alpha_est,type="l",lwd=3,ylab="alpha estimation",xlab="iterasi"),plot(beta_est,type="l",lwd=3,ylab="beta estimation",xlab="iterasi")))
}

2. Distribusi Weibull

pdf Weibull adalah \(f(x)\)= $ ()^{-1} e^$

#nr distribusi weibull
nr.weibull <- function(x,theta0){
  e=10^-6
  n=length(x)
  theta=matrix(nrow=2,ncol=200)
  theta[,1] = theta0
  error <- c()
  i=1
  repeat {
    alpha= theta[1,i]
    beta= theta[2,i]
    w1 = sum((x/alpha)^beta)
    w2 = sum((x/alpha)^beta*log(x/alpha))
    w3 = sum((x/alpha)^beta*log(x/alpha)^2)
    gradien = rbind(-n*beta/alpha+beta*w1/alpha,
                    n/beta-n*log(alpha)+sum(log(x))-w2)
    hess = matrix(c(n*beta/alpha^2-beta*(beta+1)*w1/alpha^2,
                    -n/alpha+w1/alpha-beta*w2/alpha,
                    -n/alpha+w1/alpha-beta*w2/alpha,
                    -n/beta^2-w3),ncol=2,byrow = TRUE)
    theta[,i+1]= theta[,i]- solve(hess,gradien)
    error[i+1]=abs(sum(theta[,i+1]-theta[,i]))
    if (error[i+1]<e){break}
    i=i+1
  }
  par(mfrow=c(2,1))
  return(list(cbind(t(theta[,1:length(error)]),error),alpha=alpha,beta=beta,plot(theta[1,],pch=19,ylab="alpha"),plot(theta[2,],pch=19,ylab="beta")))
}

3. Distribusi Eksponensial

#estimasi parameter eksponensial
nr.exp <- function(x,lambda0){
n <- length(x)
lambda_est=c();error_exp=c()
lambda_est[1] <- lambda0
e=10^-6
i=1
repeat{
  lambda <- lambda_est[i]
  lambda_est[i+1]<-lambda_est[i]+(lambda^2/n)*(n/lambda-sum(x))
  error_exp[i+1]<-abs(lambda_est[i+1]-lambda_est[i])
  if(error_exp[i+1]<e){break}
  i=i+1
}
return(list(iterasi=cbind(lambda_est,error_exp),parameter=lambda))

}
nr.normal <- function(x,paramawal){
e <- 10^-6
param.est <- matrix(ncol = 100, nrow = 2)
param.est[,1] <- paramawal
n <- length(x)
i=1
error <- c()
repeat
{
  miu <- param.est[1,i]
  sigma2 <- param.est[2,i]
  #membuat vektor gradien
  grad1 <- (1/sigma2)*(sum(x-miu))
  grad2 <- (-n/(2*sigma2))+(1/(2*sigma2^2)*sum((x-miu)^2))
  gradien <- as.matrix(c(grad1,grad2),nrow=2,ncol=1)
  # membuat matriks hessian
  h11 <- (-n/sigma2)
  h21 <- -(sum(x-miu))/sigma2^2
  h22 <- -(n/(2*sigma2^2))-((1/sigma2^3)*sum((x-miu)^2))
  hessian <- matrix(c(h11,h21,h21,h22), ncol = 2, nrow = 2, byrow = TRUE)
  param <- as.matrix(param.est[,i])-(solve(hessian)%*%gradien)
  param.est[,i+1] <- param
  
  error[i+1] <- (abs(sum(param.est[,i+1]-param.est[,i])))
  if (error[i+1]<e) {break}
  i=i+1
}
return(list(hasil=cbind(t(param.est[,1:length(error)]),error),estimasi=param,plot(param.est[1,],pch=19,ylab="miu"),plot(param.est[2,],pch=19,ylab="sigma2")))
}
x=c(17.23,16.18,13.94,23.23,27.08,18.53,21.22,8.23,7.42,4.87,8.55,7.60)
theta0=c(3,2)
nr.gamma(x,theta0)

## [[1]]
##      alpha_est beta_est        error
## [1,]  3.000000 2.000000           NA
## [2,]  3.767766 4.835556 7.677665e-01
## [3,]  4.028942 3.850203 2.611756e-01
## [4,]  4.049113 3.600614 2.017059e-02
## [5,]  4.049218 3.582678 1.052555e-04
## [6,]  4.049218 3.582585 2.836751e-09
## 
## $alpha
## [1] 4.049218
## 
## $beta
## [1] 3.582585
## 
## [[4]]
## NULL
## 
## [[5]]
## NULL

Estimasi parameter Gamma menunjukkan bahwa \(\alpha\) kovergen pada nilai 4.0492 dan \(\beta\) konvergen pada nilai 3.582

nr.weibull(x,theta0)

## [[1]]
##                                   error
##   [1,]  3.000000 2.0000000           NA
##   [2,]  5.592169 0.8358360 1.428005e+00
##   [3,]  7.990872 0.7087583 2.271626e+00
##   [4,] 10.103261 0.9993627 2.402993e+00
##   [5,] 12.039171 1.3075135 2.244060e+00
##   [6,] 13.578925 1.6112946 1.843535e+00
##   [7,] 14.694713 1.8623021 1.366795e+00
##   [8,] 15.446898 2.0385506 9.284331e-01
##   [9,] 15.907922 2.1492765 5.717499e-01
##  [10,] 16.170134 2.2118908 3.248264e-01
##  [11,] 16.308306 2.2454076 1.716896e-01
##  [12,] 16.379583 2.2622447 8.811380e-02
##  [13,] 16.414189 2.2708181 4.317963e-02
##  [14,] 16.431957 2.2748397 2.178893e-02
##  [15,] 16.439901 2.2769761 1.008057e-02
##  [16,] 16.444444 2.2778517 5.418576e-03
##  [17,] 16.446012 2.2784303 2.146216e-03
##  [18,] 16.447373 2.2785668 1.497573e-03
##  [19,] 16.447468 2.2787687 2.971429e-04
##  [20,] 16.448046 2.2787389 5.480772e-04
##  [21,] 16.447826 2.2788446 1.144492e-04
##  [22,] 16.448190 2.2787828 3.023859e-04
##  [23,] 16.447923 2.2788593 1.902353e-04
##  [24,] 16.448211 2.2787962 2.248348e-04
##  [25,] 16.447958 2.2788599 1.894363e-04
##  [26,] 16.448205 2.2788020 1.890759e-04
##  [27,] 16.447977 2.2788575 1.727524e-04
##  [28,] 16.448193 2.2788057 1.651407e-04
##  [29,] 16.447990 2.2788547 1.542143e-04
##  [30,] 16.448182 2.2788087 1.457556e-04
##  [31,] 16.448002 2.2788520 1.369062e-04
##  [32,] 16.448172 2.2788113 1.290085e-04
##  [33,] 16.448012 2.2788497 1.213635e-04
##  [34,] 16.448162 2.2788136 1.142712e-04
##  [35,] 16.448021 2.2788475 1.075438e-04
##  [36,] 16.448154 2.2788156 1.012375e-04
##  [37,] 16.448028 2.2788457 9.528799e-05
##  [38,] 16.448146 2.2788173 8.969521e-05
##  [39,] 16.448035 2.2788440 8.442654e-05
##  [40,] 16.448140 2.2788189 7.946999e-05
##  [41,] 16.448042 2.2788425 7.480259e-05
##  [42,] 16.448134 2.2788203 7.041071e-05
##  [43,] 16.448047 2.2788412 6.627557e-05
##  [44,] 16.448129 2.2788215 6.238421e-05
##  [45,] 16.448052 2.2788401 5.872054e-05
##  [46,] 16.448125 2.2788226 5.527272e-05
##  [47,] 16.448056 2.2788390 5.202673e-05
##  [48,] 16.448121 2.2788236 4.897190e-05
##  [49,] 16.448060 2.2788381 4.609597e-05
##  [50,] 16.448117 2.2788244 4.338935e-05
##  [51,] 16.448063 2.2788373 4.084129e-05
##  [52,] 16.448114 2.2788252 3.844319e-05
##  [53,] 16.448066 2.2788366 3.618561e-05
##  [54,] 16.448111 2.2788258 3.406086e-05
##  [55,] 16.448069 2.2788360 3.206065e-05
##  [56,] 16.448109 2.2788264 3.017810e-05
##  [57,] 16.448071 2.2788354 2.840591e-05
##  [58,] 16.448106 2.2788270 2.673795e-05
##  [59,] 16.448073 2.2788349 2.516779e-05
##  [60,] 16.448104 2.2788274 2.368996e-05
##  [61,] 16.448075 2.2788345 2.229880e-05
##  [62,] 16.448103 2.2788278 2.098943e-05
##  [63,] 16.448077 2.2788341 1.975686e-05
##  [64,] 16.448101 2.2788282 1.859674e-05
##  [65,] 16.448078 2.2788337 1.750468e-05
##  [66,] 16.448100 2.2788285 1.647681e-05
##  [67,] 16.448079 2.2788334 1.550924e-05
##  [68,] 16.448099 2.2788288 1.459854e-05
##  [69,] 16.448080 2.2788332 1.374127e-05
##  [70,] 16.448098 2.2788291 1.293438e-05
##  [71,] 16.448081 2.2788329 1.217484e-05
##  [72,] 16.448097 2.2788293 1.145993e-05
##  [73,] 16.448082 2.2788327 1.078697e-05
##  [74,] 16.448096 2.2788295 1.015356e-05
##  [75,] 16.448083 2.2788325 9.557317e-06
##  [76,] 16.448095 2.2788297 8.996107e-06
##  [77,] 16.448084 2.2788324 8.467834e-06
##  [78,] 16.448094 2.2788298 7.970598e-06
##  [79,] 16.448084 2.2788322 7.502547e-06
##  [80,] 16.448094 2.2788300 7.061992e-06
##  [81,] 16.448085 2.2788321 6.647296e-06
##  [82,] 16.448093 2.2788301 6.256962e-06
##  [83,] 16.448086 2.2788320 5.889540e-06
##  [84,] 16.448093 2.2788302 5.543701e-06
##  [85,] 16.448086 2.2788319 5.218164e-06
##  [86,] 16.448092 2.2788303 4.911748e-06
##  [87,] 16.448086 2.2788318 4.623321e-06
##  [88,] 16.448092 2.2788304 4.351835e-06
##  [89,] 16.448087 2.2788317 4.096287e-06
##  [90,] 16.448092 2.2788305 3.855749e-06
##  [91,] 16.448087 2.2788316 3.629332e-06
##  [92,] 16.448091 2.2788305 3.416214e-06
##  [93,] 16.448087 2.2788316 3.215608e-06
##  [94,] 16.448091 2.2788306 3.026783e-06
##  [95,] 16.448087 2.2788315 2.849045e-06
##  [96,] 16.448091 2.2788307 2.681746e-06
##  [97,] 16.448088 2.2788315 2.524269e-06
##  [98,] 16.448091 2.2788307 2.376041e-06
##  [99,] 16.448088 2.2788314 2.236516e-06
## [100,] 16.448091 2.2788307 2.105185e-06
## [101,] 16.448088 2.2788314 1.981565e-06
## [102,] 16.448090 2.2788308 1.865205e-06
## [103,] 16.448088 2.2788313 1.755677e-06
## [104,] 16.448090 2.2788308 1.652582e-06
## [105,] 16.448088 2.2788313 1.555539e-06
## [106,] 16.448090 2.2788308 1.464196e-06
## [107,] 16.448088 2.2788313 1.378216e-06
## [108,] 16.448090 2.2788309 1.297286e-06
## [109,] 16.448088 2.2788313 1.221107e-06
## [110,] 16.448090 2.2788309 1.149402e-06
## [111,] 16.448089 2.2788312 1.081907e-06
## [112,] 16.448090 2.2788309 1.018376e-06
## [113,] 16.448089 2.2788312 9.585755e-07
## 
## $alpha
## [1] 16.44809
## 
## $beta
## [1] 2.278831
## 
## [[4]]
## NULL
## 
## [[5]]
## NULL

Estimasi parameter Distribusi Weibull menggunakan NR menunjukkan bahwa \(\alpha\) konvergen pada nilai 16.44 dan \(\beta\) konvergen pada nilai 2.278.

nr.normal(x,theta0)

## $hasil
##                                   error
##  [1,]  3.00000  2.00000000           NA
##  [2,] 24.77033  0.21604969 1.998638e+01
##  [3,]  8.28103  0.08500026 1.662035e+01
##  [4,] 13.93257  0.09283859 5.659376e+00
##  [5,] 14.22244  0.13880146 3.358348e-01
##  [6,] 14.36541  0.20778179 2.119547e-01
##  [7,] 14.43653  0.31095102 1.742862e-01
##  [8,] 14.47194  0.46489319 1.893559e-01
##  [9,] 14.48956  0.69394776 2.466689e-01
## [10,] 14.49830  1.03339261 3.481843e-01
## [11,] 14.50262  1.53345552 5.043817e-01
## [12,] 14.50474  2.26374762 7.324130e-01
## [13,] 14.50577  3.31681280 1.054097e+00
## [14,] 14.50626  4.80784228 1.491523e+00
## [15,] 14.50649  6.86531771 2.057706e+00
## [16,] 14.50660  9.60579758 2.740584e+00
## [17,] 14.50664 13.08902702 3.483273e+00
## [18,] 14.50666 17.26203716 4.173027e+00
## [19,] 14.50666 21.92129879 4.659267e+00
## [20,] 14.50667 26.73141732 4.810120e+00
## [21,] 14.50667 31.31185570 4.580439e+00
## [22,] 14.50667 35.35001452 4.038159e+00
## [23,] 14.50667 38.67520019 3.325186e+00
## [24,] 14.50667 41.26277857 2.587578e+00
## [25,] 14.50667 43.18917002 1.926391e+00
## [26,] 14.50667 44.57656927 1.387399e+00
## [27,] 14.50667 45.55210640 9.755371e-01
## [28,] 14.50667 46.22653699 6.744306e-01
## [29,] 14.50667 46.68736290 4.608259e-01
## [30,] 14.50667 46.99971866 3.123558e-01
## [31,] 14.50667 47.21028821 2.105696e-01
## [32,] 14.50667 47.35171927 1.414311e-01
## [33,] 14.50667 47.44647836 9.475909e-02
## [34,] 14.50667 47.50986207 6.338371e-02
## [35,] 14.50667 47.55221204 4.234997e-02
## [36,] 14.50667 47.58048732 2.827528e-02
## [37,] 14.50667 47.59935620 1.886888e-02
## [38,] 14.50667 47.61194377 1.258757e-02
## [39,] 14.50667 47.62033918 8.395411e-03
## [40,] 14.50667 47.62593777 5.598586e-03
## [41,] 14.50667 47.62967089 3.733122e-03
## [42,] 14.50667 47.63215996 2.489073e-03
## [43,] 14.50667 47.63381949 1.659527e-03
## [44,] 14.50667 47.63492590 1.106415e-03
## [45,] 14.50667 47.63566354 7.376388e-04
## [46,] 14.50667 47.63615531 4.917719e-04
## [47,] 14.50667 47.63648317 3.278536e-04
## [48,] 14.50667 47.63670174 2.185716e-04
## [49,] 14.50667 47.63684746 1.457155e-04
## [50,] 14.50667 47.63694460 9.714415e-05
## [51,] 14.50667 47.63700936 6.476299e-05
## [52,] 14.50667 47.63705254 4.317542e-05
## [53,] 14.50667 47.63708132 2.878366e-05
## [54,] 14.50667 47.63710051 1.918912e-05
## [55,] 14.50667 47.63711330 1.279276e-05
## [56,] 14.50667 47.63712183 8.528509e-06
## [57,] 14.50667 47.63712752 5.685675e-06
## [58,] 14.50667 47.63713131 3.790451e-06
## [59,] 14.50667 47.63713383 2.526967e-06
## [60,] 14.50667 47.63713552 1.684645e-06
## [61,] 14.50667 47.63713664 1.123097e-06
## [62,] 14.50667 47.63713739 7.487312e-07
## 
## $estimasi
##          [,1]
## [1,] 14.50667
## [2,] 47.63714
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL
lambda0= seq(0.01, 0.5, length = 100)
exp=c()
lambda_opt=c()
for (i in 1:length(lambda0)){
  lambda_opt[i] = nr.exp(x,lambda0[1])$parameter
}
lambda_opt
##   [1] 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382
##   [7] 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382
##  [13] 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382
##  [19] 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382
##  [25] 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382
##  [31] 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382
##  [37] 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382
##  [43] 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382
##  [49] 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382
##  [55] 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382
##  [61] 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382
##  [67] 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382
##  [73] 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382
##  [79] 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382
##  [85] 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382
##  [91] 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382 0.06893382
##  [97] 0.06893382 0.06893382 0.06893382 0.06893382
#plot mendapatkan lambda optimum
n=length(x)
N=length(lambda0)
ln=c(1:N)
for (i in 1:N)
{ln[i]=prod(dexp(x,lambda0[i]))}
plot(lambda0,ln,main="parameter optimum dist. Eksponensial",xlab="lambda",ylab="ln Likelihood")

cat("nilai lambda optimum ",lambda0[which.max(ln)])
## nilai lambda optimum  0.06939394

Optimasi parameter menggunakan optim

Untuk membuktikan kebenaran hasil NR maka digunakan function optim

#optimasi untuk parameter gamma
gamma.opt <- function(theta,x){
  pdf=matrix(nrow=length(x),ncol=2)
  alpha <- theta[1]
  beta <- theta[2]
  pdf[,1] <- dgamma(x,shape=alpha,scale=beta)
  pdf[,2] <- log10(pdf[,1])
  sumlog <- -sum(pdf[,2])
  fn <- function(theta){
    a <- theta[1]
    b <- theta[2]
    fy <- dgamma(x,shape=a,scale=b)
    logfy <- log10(fy)
    sumlog1 <- -sum(logfy)
  }
  theta.opt <- optim(theta,fn)
  theta_opt <- theta.opt$par
  pdff <- matrix(nrow=length(x),ncol=2)
  pdff[,1] <- dgamma(x,shape=theta_opt[1],scale=theta_opt[2])
  pdff[,2] <- log10(pdff[,1])
  sumlogg <- -sum(pdf[,2])
  pdff1 <- cbind(pdf,pdff)
  hasil <- list(pdf=pdff1,sumlog=sumlogg,theta=theta_opt)
  print(hasil)
}
gamma.opt(theta0,x)
## $pdf
##               [,1]       [,2]       [,3]      [,4]
##  [1,] 3.365144e-03 -2.4729964 0.04283366 -1.368215
##  [2,] 5.016429e-03 -2.2996053 0.04740441 -1.324181
##  [3,] 1.141230e-02 -1.9426270 0.05624107 -1.249946
##  [4,] 3.045426e-04 -3.5163520 0.01995456 -1.699958
##  [5,] 6.037123e-05 -4.2191699 0.01087308 -1.963647
##  [6,] 2.031854e-03 -2.6921075 0.03719659 -1.429497
##  [7,] 6.942347e-04 -3.1584937 0.02653892 -1.576117
##  [8,] 6.911271e-02 -1.1604421 0.05551891 -1.255559
##  [9,] 8.422777e-02 -1.0745447 0.05074983 -1.294565
## [10,] 1.298467e-01 -0.8865691 0.02863755 -1.543064
## [11,] 6.356285e-02 -1.1967967 0.05703761 -1.243839
## [12,] 8.075849e-02 -1.0928118 0.05192239 -1.284645
## 
## $sumlog
## [1] 25.71252
## 
## $theta
## [1] 4.049432 3.582059
#optimasi untuk parameter weibul
weibull.opt <- function(theta,x){
  pdf=matrix(nrow=length(x),ncol=2)
  alpha <- theta[2]
  beta <- theta[1]
  pdf[,1] <- dweibull(x,shape=alpha,scale=beta)
  pdf[,2] <- log10(pdf[,1])
  sumlog <- -sum(pdf[,2])
  fn <- function(theta){
    a <- theta[2]
    b <- theta[1]
    fy <- dweibull(x,shape=a,scale=b)
    logfy <- log10(fy)
    sumlog1 <- -sum(logfy)
  }
  theta.opt <- optim(theta,fn)
  theta_opt <- theta.opt$par
  pdff <- matrix(nrow=length(x),ncol=2)
  pdff[,1] <- dweibull(x,shape=theta_opt[2],scale=theta_opt[1])
  pdff[,2] <- log10(pdff[,1])
  sumlogg <- -sum(pdf[,2])
  pdff1 <- cbind(pdf,pdff)
  hasil <- list(pdf=pdff1,sumlog=sumlogg,theta=theta_opt)
  print(hasil)
}
weibull.opt(theta0,x)
## Warning in dweibull(x, shape = a, scale = b): NaNs produced

## Warning in dweibull(x, shape = a, scale = b): NaNs produced
## $pdf
##               [,1]       [,2]       [,3]      [,4]
##  [1,] 1.809206e-14 -13.742512 0.04839544 -1.315196
##  [2,] 8.375077e-13 -12.077011 0.05179519 -1.285711
##  [3,] 1.300098e-09  -8.886024 0.05647794 -1.248121
##  [4,] 4.708496e-26 -25.327118 0.02398236 -1.620108
##  [5,] 2.470660e-35 -34.607187 0.01164430 -1.933886
##  [6,] 1.111226e-16 -15.954198 0.04347197 -1.361791
##  [7,] 8.808299e-22 -21.055108 0.03216158 -1.492663
##  [8,] 9.856895e-04  -3.006260 0.04646439 -1.332880
##  [9,] 3.634528e-03  -2.439552 0.04249852 -1.371626
## [10,] 7.759948e-02  -1.110141 0.02741423 -1.562024
## [11,] 5.638930e-04  -3.248803 0.04788352 -1.319814
## [12,] 2.756739e-03  -2.559604 0.04342347 -1.362275
## 
## $sumlog
## [1] 144.0135
## 
## $theta
## [1] 16.452914  2.279682
#menggunakan optimasi untuk ekponensial
exp.opt <- function(lambda,x){
  pdf=matrix(nrow=length(x),ncol=2)
  pdf[,1] <- (lambda)*exp(-x*lambda)
  pdf[,2] <- log10(pdf[,1])
  sumlog <- -sum(pdf[,2])
  fn <- function(xx){
    y <-xx[1]
    fy <- (y)*exp(-x*y)
    logfy <- log10(fy)
    sumlog1 <- -sum(logfy)
  }
  lambda.opt <- optim(lambda,fn,method="BFGS")
  lambda.opt <- lambda.opt$par
  pdff <- matrix(nrow=length(x),ncol=2)
  pdff[,1] <- (lambda.opt)*exp(-x*lambda.opt)
  pdff[,2] <- log10(pdff[,1])
  sumlogg <- -sum(pdf[,2])
  pdff1 <- cbind(pdf,pdff)
  hasil <- list(pdf=pdff1,sumlog=sumlogg,lambda=lambda.opt)
  print(hasil)
}
exp.opt(0.2,x)
## Warning in fn(par, ...): NaNs produced

## Warning in fn(par, ...): NaNs produced

## Warning in fn(par, ...): NaNs produced

## Warning in fn(par, ...): NaNs produced

## Warning in fn(par, ...): NaNs produced
## $pdf
##               [,1]      [,2]       [,3]      [,4]
##  [1,] 0.0063745746 -2.195549 0.02101892 -1.677390
##  [2,] 0.0078641729 -2.104347 0.02259664 -1.645956
##  [3,] 0.0123088359 -1.909783 0.02636951 -1.578898
##  [4,] 0.0019199850 -2.716702 0.01389922 -1.857010
##  [5,] 0.0008889781 -3.051109 0.01065944 -1.972266
##  [6,] 0.0049151259 -2.308465 0.01921731 -1.716307
##  [7,] 0.0028700153 -2.542116 0.01596479 -1.796837
##  [8,] 0.0385639293 -1.413819 0.03908772 -1.407960
##  [9,] 0.0453457913 -1.343463 0.04133225 -1.383711
## [10,] 0.0755139469 -1.121973 0.04927510 -1.307372
## [11,] 0.0361731585 -1.441614 0.03823496 -1.417539
## [12,] 0.0437423774 -1.359098 0.04082258 -1.389100
## 
## $sumlog
## [1] 23.50804
## 
## $lambda
## [1] 0.06893172
norm.opt <- function(x,theta){
  pdf=matrix(nrow=length(x),ncol=2)
  miu <- theta[1]
  sigma2 <- theta[2]
  pdf[,1] <- dnorm(x,miu,sigma2)
  pdf[,2] <- log10(pdf[,1])
  sumlog <- -sum(pdf[,2])
  fn <- function(theta){
    miu2 <- theta[1]
    sigma22 <- theta[2]
    fy <- dnorm(x,miu2,sigma22)
    logfy <- log10(fy)
    sumlog1 <- -sum(logfy)
  }
  theta.opt <- optim(theta,fn)
  theta_opt <- theta.opt$par
  pdff <- matrix(nrow=length(x),ncol=2)
  pdff[,1] <- dnorm(x,theta_opt[1],theta_opt[2])
  pdff[,2] <- log10(pdff[,1])
  sumlogg <- -sum(pdf[,2])
  pdff1 <- cbind(pdf,pdff)
  hasil <- list(pdf=pdff1,sumlog=sumlogg,theta=theta_opt)
  return(hasil)
}
norm.opt(x,theta0)
## $pdf
##               [,1]        [,2]       [,3]      [,4]
##  [1,] 2.028553e-12 -11.6928136 0.05348032 -1.271806
##  [2,] 7.406073e-11 -10.1304120 0.05613839 -1.250740
##  [3,] 6.348033e-08  -7.1973608 0.05762224 -1.239410
##  [4,] 1.210168e-23 -22.9171545 0.02599559 -1.585100
##  [5,] 6.635423e-33 -32.1781314 0.01098794 -1.959084
##  [6,] 1.610417e-14 -13.7930617 0.04877300 -1.311821
##  [7,] 1.898342e-19 -18.7216255 0.03601054 -1.443570
##  [8,] 6.530980e-03  -2.1850216 0.03823328 -1.417558
##  [9,] 1.735047e-02  -1.7606888 0.03412590 -1.466916
## [10,] 1.288382e-01  -0.8899555 0.02180838 -1.661377
## [11,] 4.243367e-03  -2.3722894 0.03983746 -1.399708
## [12,] 1.416352e-02  -1.8488288 0.03504044 -1.455430
## 
## $sumlog
## [1] 125.6873
## 
## $theta
## [1] 14.505466  6.900199
cat("nilai varians adalah ",norm.opt(x,theta0)$theta[2]^2)
## nilai varians adalah  47.61275