1

如此处所述,这些madness软件包非常适合 R 中的 autodiff

我现在想计算一个导数 wrt y 的导数 wrt x。

$\frac{\partial}{\partial x}\frac{\partial}{\partial y}xy$

如何使用madness

更新:实际上我猜这是因素..也许只需将两个导数相乘就可以了吗?也许只有当 x 是 y 的函数时,这才会变得困难。

4

1 回答 1

1

这是使用numderiv函数 in的一种方法madness

library(madness)

dxdy <- function(x, y, f) {
  dy <- function(x, y) {
    dvdx(f(x, y))
  }
  
  numderiv(dy, x, y = y)
}

x <- matrix(1:3, nrow = 1)
y <- matrix(3:1, ncol = 1)
# identity matrix, as expected
dvdx(dxdy(madness(x), madness(y), function(x, y) x%*%y))
#>      [,1] [,2] [,3]
#> [1,]    1    0    0
#> [2,]    0    1    0
#> [3,]    0    0    1
x <- matrix(2, ncol = 1)
y <- matrix(3, ncol = 1)
dvdx(dxdy(madness(x), madness(y), function(x, y) y^x))
#>          [,1]
#> [1,] 9.591674
# compare to analytical solution
y^(x-1)*(x*log(y) + 1)
#>          [,1]
#> [1,] 9.591674
x <- matrix(1:3, ncol = 1)
y <- matrix(3:1, ncol = 1)
dvdx(dxdy(madness(x), madness(y), function(x, y) sum(y^x)))
#>          [,1]     [,2] [,3]
#> [1,] 2.098612 0.000000    0
#> [2,] 0.000000 4.772589    0
#> [3,] 0.000000 0.000000    1
# compare to analytical solution
y^(x-1)*(x*log(y) + 1)
#>          [,1]
#> [1,] 2.098612
#> [2,] 4.772589
#> [3,] 1.000000
x <- matrix(1:3, ncol = 1)
y <- matrix(3:1, ncol = 1)
dvdx(dxdy(madness(x), madness(y), function(x, y) sum(sin(x*y))))
#>           [,1]     [,2]      [,3]
#> [1,] -1.413352 0.000000  0.000000
#> [2,]  0.000000 2.373566  0.000000
#> [3,]  0.000000 0.000000 -1.413353
# compare to analytical solution
cos(x*y) - x*y*sin(x*y)
#>           [,1]
#> [1,] -1.413353
#> [2,]  2.373566
#> [3,] -1.413353

请注意,madness导数是每个x(列)相对于y每个xy(行)的部分的部分。

于 2021-11-10T14:53:40.610 回答