先將公式擺好
options(digits = 5)
rm(list=ls())
MLE <- function(gamma,psi){
# equation 12
dfdg <-
(512*gamma^3 + 1280*gamma^2 + 1600*gamma*psi - 800*psi - 3200) /
((16*gamma^2 + 10*psi + 40)^2)
# equation 15
dfdp <-
(-480*gamma^2 + 800*gamma + 100*psi - 1600) /
((16*gamma^2 + 10*psi + 40)^2)
# equation 19
df2dg2 <-
(-8192*gamma^4 - 40960*gamma^3 + 61440*gamma^2 + 307200*gamma - 61440*gamma^2*psi + 76800*gamma*psi +16000*psi^2 + 64000*psi) /
((16*gamma^2 + 10*psi + 40)^3)
# equation 22
df2dp2 <-
(11200*gamma^2 - 16000*gamma - 1000*psi + 36000) /
((16*gamma^2 + 10*psi + 40)^3)
# equation 25
df2dgp <-
(15360*gamma^3 - 38400*gamma^2 + 64000*gamma - 16000*gamma*psi + 8000*psi +32000) /
((16*gamma^2 + 10*psi + 40)^3)
# 目標函數
FML <- log(16*gamma^2 + 10*psi + 40) +
((16*gamma^2 + 10*psi + 40)^-1) * (80*gamma^2-80*gamma + 10*psi + 240) -
log(175) - 2
# matirx
matgp <- matrix(c(gamma,psi),nrow = 2)
matf1 <- matrix(c(df2dg2,df2dgp,df2dgp,df2dp2),nrow=2)
matf2 <- matrix(c(dfdg,dfdp),nrow=2)
new.mat<- matgp - solve(matf1) %*% matf2
sol<- c(new.mat,FML)
names(sol) <- c("gamma","psi","FML")
sol
}
# 以(0.5,5)當起始值跟課本 equation 28 一致
MLE(0.5,5)
## gamma psi FML
## 0.61015 7.74715 0.25085
寫個迴圈來iteration吧
Iteration <- function(times,ini_gamma,ini_psi){
for (i in 1L:times){
tab <- c(i,MLE(ini_gamma,ini_psi))
names(tab) <- c("iterations","gamma","psi","FML")
mat <- MLE(ini_gamma,ini_psi)[1:2]
ini_gamma <- mat[1]
ini_psi <- mat[2]
print(tab)
}}
Iteration(15L,0.5,5)
## iterations gamma psi FML
## 1.00000 0.61015 7.74715 0.25085
## iterations gamma psi FML
## 2.000000 0.622959 10.530680 0.068817
## iterations gamma psi FML
## 3.00000 0.62483 12.31961 0.01090
## iterations gamma psi FML
## 4.0000e+00 6.2500e-01 1.2841e+01 5.2636e-04
## iterations gamma psi FML
## 5.0000e+00 6.2500e-01 1.2875e+01 1.9156e-06
## iterations gamma psi FML
## 6.0000e+00 6.2500e-01 1.2875e+01 2.9089e-11
## iterations gamma psi FML
## 7.000 0.625 12.875 0.000
## iterations gamma psi FML
## 8.000 0.625 12.875 0.000
## iterations gamma psi FML
## 9.000 0.625 12.875 0.000
## iterations gamma psi FML
## 10.000 0.625 12.875 0.000
## iterations gamma psi FML
## 11.000 0.625 12.875 0.000
## iterations gamma psi FML
## 12.000 0.625 12.875 0.000
## iterations gamma psi FML
## 13.000 0.625 12.875 0.000
## iterations gamma psi FML
## 14.000 0.625 12.875 0.000
## iterations gamma psi FML
## 15.000 0.625 12.875 0.000
視覺化迭代過程
#devtools::install_github("dgrtwo/gganimate")
# load required packages
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(readr)
library(ggplot2)
library(gganimate)
library(scales)
##
## Attaching package: 'scales'
## The following object is masked from 'package:readr':
##
## col_factor
library(gapminder)
ini_gamma <- 0.5
ini_psi <- 5
tab.all <- c()
for (i in 1L:15L){
tab <- c(i,MLE(ini_gamma,ini_psi))
names(tab) <- c("iterations","gamma","psi","FML")
mat <- MLE(ini_gamma,ini_psi)[1:2]
ini_gamma <- mat[1]
ini_psi <- mat[2]
tab.all <- c(tab.all,tab)
}
tab.all<- tab.all %>% matrix(.,ncol=4,byrow = T) %>% data.frame()
names(tab.all) <- c("iterations","gamma","psi","FML")
tab.all$iterations<- as.integer(tab.all$iterations)
p1 <- ggplot(tab.all, aes(gamma, psi, colour = iterations)) +
geom_point(alpha = 0.7) +
scale_size(range = c(0, 3)) +
theme(legend.position = 'none') +
labs(title = 'iterations: {frame_time}', x = 'gamma', y = 'psi') +
transition_time(iterations) +
ease_aes('linear')
animate(p1, 100, 10)
