daily pastebin goal
46%
SHARE
TWEET

Untitled

a guest Mar 22nd, 2019 77 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. # 单对初始值
  2. MultiD.Newton.Raphson <- function(expr, namevec, x0,dx = NULL, max.iter = 10000, c = 1e-09){
  3.   # 临界条件
  4.   if(length(namevec) != length(x0)){stop("namevec must have the same length with x0")}
  5.   judge <- function(x1,x0,value){
  6.     ## a problem:input initial value vector and eval derivative function##
  7.     value1 <- do.call(dx,as.list(x1))
  8.     dx1 <- attributes(value1)$gradient
  9.     dx0 <- attributes(value)$gradient
  10.     if(abs(value1-value)<c){return(TRUE)} #f(x1)-f(x0)
  11.     if(max(abs(x1-x0))<c){return(TRUE)} # x1-x0
  12.     if(max(abs(dx1))<c){return(TRUE)} #f'(x1) - f'(x0)
  13.     return(FALSE)
  14.   }
  15.   # 判断海塞矩阵是否正定、负定还是不定,以确定极小/大值
  16.   is.positive.definite <- function(x){
  17.     eigens <- eigen(x)$values
  18.     if(length(eigens>0) == length(eigens)){
  19.       return('min')
  20.     }
  21.     if(length(eigens<0) == length(eigens)){
  22.       return('max')
  23.     }else{
  24.       return('不是极值')
  25.     }
  26.   }
  27.  
  28.   if(is.null(dx)){
  29.     dx <- deriv(expr, namevec,hessian = TRUE, function.arg = TRUE)
  30.     for (k in 1:max.iter) {
  31.       value <- do.call(dx, as.list(x0)) # pass a vector to  function with more than one arguments
  32.       A <- rbind(attributes(value)$hessian[,,1:length(namevec)]) # Hessian matrix
  33.       b <- attributes(value)$gradient #b
  34.       stepsize <- solve(A,t(b))
  35.       x1 <- x0 - stepsize
  36.       if(judge(x1,x0,value)){ # 临界条件
  37.         #返回迭代次数、向量X的值、海塞矩阵、极值、是极大值还是极小值
  38.         return(list(it.num = k, x = x1, hessian = A, pv = as.numeric(value), pole = is.positive.definite(A)))
  39.       }
  40.       x0 <- x1
  41.     }
  42.   }else{
  43.     # if dx exists, it must be a list that has hessian matrix and derivae
  44.     xlist <- setNames(as.list(x0), namevec)
  45.     fx <- eval(expr, envir = xlist)
  46.     A <- do.call(dx, as.list(x0))$hessian
  47.     b <- do.call(dx, as.list(x0))$gradient
  48.     stepsize <- solve(A,t(b))
  49.     x1 <- x0 - stepsize
  50.     if(judge(x1,x0,value)){ # 临界条件
  51.       #返回迭代次数、向量X的值、海塞矩阵、极值、是极大值还是极小值
  52.       return(list(it.num = k, x = x1, hessian = A, pv = as.numeric(value), pole = is.positive.definite(A)))
  53.     }
  54.     x0 <- x1
  55.   }
  56.   return(list(it.num = max.iter, x = x0, hessian = A, pv = as.numeric(value), pole = is.positive.definite(A)))
  57. }
  58.  
  59. # 适用于多个初始值的情况,要求输入为X Y列表
  60. MultiD.Newton.Raphson.initvalues <- function(expr, namevec, x0, aim, dx=NULL, max.iter = 10000, c = 1e-09){
  61.   # x0 must be a  list
  62.   pv <- list()
  63.   for (i in 1:length(x0)) {
  64.     pv[[i]] <- MultiD.Newton.Raphson(expr,namevec, x0[[i]],dx,max.iter,c)
  65.   }
  66.   # aim ==1 max, aim == 0, min
  67.   if(aim == 1){
  68.     valus <- sapply(pv, function(x){return(x$pv)})
  69.     return(pv[[which(valus == max(valus))]])
  70.   }
  71.   if(aim == 0){
  72.     valus <- sapply(pv, function(x){return(x$pv)})
  73.     return(pv[[which(valus == min(valus))]])
  74.   }
  75. }
  76.  
  77. MultiD.Newton.Raphson.initvalues(expression(x^2-x*y+y^2+exp(y)),
  78.                                  namevec = c('x','y'),
  79.                                  x0 = list(2:3,3:4),
  80.                                  aim = 0)
  81.  
  82. #$it.num
  83. #[1] 7
  84.  
  85. #$x
  86. #[,1]
  87. #x -0.2162814
  88. #y -0.4325628
  89.  
  90. #$hessian
  91. #x        y
  92. #x  2 -1.00000
  93. #y -1  2.64885
  94.  
  95. #$pv
  96. #[1] 0.789177
  97.  
  98. #$pole
  99. #[1] "min"
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top