• API
• FAQ
• Tools
• Archive
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))
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
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
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.

Top