计算多项式的导数和曲率

Question

我有一个数据框，记录了一段时间内一群人的皮肤温度。我愿意：

在 Time 上对每个 id 的 SkinTemp 拟合二次多项式；
计算曲率。

这似乎比它应该的要难得多。

我在中问了第一部分，但是我无法继续计算导数和曲率。

df <- data.frame(Time = seq(65),
                 SkinTemp = rnorm(65, 37, 0.5),
                 id = rep(1:10, c(5,4,10,6,7,8,9,8,4,4)))

#Predict data points for each quadratic
fitted_models = df %>% group_by(id) %>% do(model =
                lm(SkinTemp ~ Time+I(Time^2), data = .))

现在我需要计算曲率 k = y''/(1 + y' ^ 2) ^ (3 / 2)，其中 y' 和 y'' 是 y 相对于 x 的一阶和二阶导数。

我想我可以要求 predict 函数通过传递例如 deriv = 2 来给我导数，但它似乎不起作用。

predQ <- lapply(unique(df$id),
                function(x) predict(deriv = 2,fitted_models$model[[x]]))

所以我修改了这个功能，它似乎工作正常但是没有这个任务的内置功能吗？

deriv <- function(x, y) diff(y) / diff(x)
middle_pts <- function(x) x[-1] - diff(x) / 2
second_d <- lapply(unique(df$id),
                   function(x) deriv(middle_pts(df$Time[df["id"]==x]), deriv(df$Time[df["id"]==x], df$SkinTemp[df["id"]==x])))

Answer 1

你只是被多项式的导数计算阻碍了。使用中定义的函数 g 怎么样？

例如，您的第一个模型具有多项式系数：

pc <- coef(fitted_models[[2]][[1]])
#(Intercept)        Time   I(Time^2) 
#38.36702120 -0.61025716  0.04703084

假设您只想评估观察位置的导数和曲率：

x <- with(df, Time[id == 1])
#[1] 1 2 3 4 5

然后就可以一步步做解析计算了：

## 1st derivative
d1 <- g(x, pc, 1)
#[1] -0.5161955 -0.4221338 -0.3280721 -0.2340104 -0.1399487

## 2nd derivative
d2 <- g(x, pc, 2)
#[1] 0.09406168 0.09406168 0.09406168 0.09406168 0.09406168

## curvature: d2 / (1 + d1 * d1) ^ (3 / 2)
d2 / (1 + d1 * d1) ^ (3 / 2)
#[1] 0.06599738 0.07355055 0.08069004 0.08683238 0.09136444

这不是比你的有限差分近似好很多吗？

注意 g 也可以计算 nderiv = 0L，即多项式本身：

g(x, pc, 0)
#[1] 37.80379 37.33463 36.95953 36.67849 36.49151

同意predict.lm:

predict.lm(fitted_models[[2]][[1]], data.frame(Time = x))
#       1        2        3        4        5 
#37.80379 37.33463 36.95953 36.67849 36.49151

函数g通过多项式系数向量pc的长度判断多项式的次数。长度为 3 的向量表示度数 = 2。它是为原始多项式设计的，而不是正交的。

所有组的计算

要对所有组进行曲率计算，我会使用 Map。

polynom_curvature <- function (x, pc) {
  d1 <- g(x, pc, 1L)
  d2 <- g(x, pc, 2L)
  d2 / (1 + d1 * d1) ^ (3 / 2)
  }

pc_lst <- lapply(fitted_models[[2]], coef)
Time_lst <- split(df$Time, df$id)
result <- Map(polynom_curvature, Time_lst, pc_lst)
str(result)
#List of 10
# $ 1 : num [1:5] 0.066 0.0736 0.0807 0.0868 0.0914
# $ 2 : num [1:4] -0.106 -0.12 -0.131 -0.135
# $ 3 : num [1:10] 0.0795 0.0897 0.0988 0.1058 0.1095 ...
# $ 4 : num [1:6] -0.098 -0.107 -0.113 -0.115 -0.112 ...
# $ 5 : num [1:7] -0.0878 -0.0923 -0.0946 -0.0944 -0.0917 ...
# $ 6 : num [1:8] 0.0752 0.0811 0.0857 0.0886 0.0895 ...
# $ 7 : num [1:9] 0.0397 0.0405 0.0411 0.0414 0.0416 ...
# $ 8 : num [1:8] 0.0178 0.018 0.0182 0.0184 0.0185 ...
# $ 9 : num [1:4] -0.151 -0.161 -0.159 -0.146
# $ 10: num [1:4] 0.1186 0.1129 0.1033 0.0917

计算多项式的导数和曲率

Calculate derivatives & curvature of a polynomial

regression

r

linear-regression

lm

polynomials

所有组的计算