【发布时间】：2021-05-12 23:36:32
【问题描述】：
总结
在我的设置中，我从给定的分隔空间中的三点轨迹开始。我将点的 x 坐标传递给一个惩罚函数 trj_points_penalty_fun，它计算轨迹的惩罚为：
-
给定三个点以及固定起点和终点的总轨迹长度
-
这些点的双变量密度值的总和。双变量以相同的定界空间为中心，具有固定的标准差和相关系数。
函数trj_points_penalty_fun 的构建方式是点的x 坐标被一个接一个地传递。然后我使用 optim() 优化函数，结果给出了惩罚最小的 x 坐标（我保持 y 坐标不变）。
问题在于，当我将点作为向量传递给惩罚函数时，如trj_vector_min_pen_fun，optim() 给我的答案与我传递时给出的答案不同一点一点。
可重现的示例
请参阅以下示例。我为篇幅道歉；我试图简化设置，但无法重现错误。
library(raster)
#> Loading required package: sp
下面的函数points_distance用于计算总轨迹的长度
points_distance <- function(x) {
x$distance <- NA
for (i in 1:(nrow(x) - 1)) {
x$distance[i] <- pointDistance(c(x[, 1][i], x[, 2][i]),
c(x[, 1][i + 1], x[, 2][i + 1]),
lonlat = FALSE
)
}
# Return output
return(x)
}
以下函数add_start_end 将固定起点和终点添加到轨迹的坐标：
add_start_end <- function(points) {
points_plus_fixed_start <- rbind(c(0,-3),points)
points_plus_fixed_start_end <- rbind(points_plus_fixed_start, c(0,3))
return (points_plus_fixed_start_end)
}
下面是一个构造函数，它创建了trj_points_penalty_fun，这是一个罚分函数，其中点数被逐一传递。
在构造函数中，双变量参数和 y 坐标作为data 传入。 x 坐标作为三个单个值传入：
MAKE_trj_points_penalty_fun <- function(data, fixed = c(FALSE, FALSE, FALSE)) {
params <- fixed
function(p) {
params[!fixed] <- p
mu <- c(0, 0)
mu[1] <- data[1]
mu[2] <- data[2]
sd1 <- data[3]
sd2 <- data[4]
rho <- data[5]
lambda <- data[6]
y_coord1 <- data[7]
y_coord2 <- data[8]
y_coord3 <- data[9]
#create vector with three points for looping.
y_coords <- c(y_coord1, y_coord2, y_coord3)
tot_density <- 0
# calculates the bivariate density for each point
for (i in seq(seq_along(params))) { # seq(1:length(params))
Q <- (params[i] - mu[1])^2 / sd1^2 + (y_coords[i] - mu[2])^2 / sd2^2 -
2 * rho * (params[i] - mu[1]) * (y_coords[i] - mu[2]) / (sd1 * sd2)
density <- 1 / (2 * pi * sd1 * sd2 * sqrt(1 - rho^2)) * exp(-Q / (2 * (1 - rho^2)))
tot_density <- tot_density + density
tot_density
}
points <- data.frame(
x = c(params[1], params[2], params[3]),
y = c(y_coord1, y_coord2, y_coord3)
)
segment <- add_start_end(points)
dist <- points_distance(segment)
trj <- sum(dist$distance, na.rm = TRUE)
# tot penalty
tot <- lambda * (trj / 10) + (1 - lambda) * (tot_density / 0.1591549)
return(tot)
}
}
这里我使用构造函数创建trj_points_penalty_fun：
trj_points_penalty_fun <- MAKE_trj_points_penalty_fun(c(0, 0, 1, 1, 0, 0.5,-2.5, 0, 2.5))
现在我们传入 x 坐标：
trj_points_penalty_fun(c(0.5, 1, 0.5))
#> [1] 0.6677013
优化功能...
optim(par = c(0, 0, 0), trj_points_penalty_fun)$par
...给出以下结果：
#> [1] 0.7146940 2.4769463 0.7148109
现在是另一种情况，其中轨迹点作为向量传递：
# Toy trajectory
toy_trj_vs <- data.frame(
x <- c(0.5, 1, 0.5),
y <- c(-2.5, 0, 2.5)
)
x_coord <- toy_trj_vs[,1]
y_coord <- toy_trj_vs[,2]
vect_f <- rep(F, length(x_coord))
以下是一个构造函数，它创建trj_vector_penalty_fun，一个惩罚函数，其中
点的坐标作为向量传递。
MAKE_trj_vector_penalty_fun <- function(data_vs, fixed_y_coord, fixed_vs =
vect_f) {
params_vs <- fixed_vs
function(p_vs) {
params_vs[!fixed_vs] <- p_vs
mu_vs <- c(0, 0)
mu_vs[1] <- data_vs[1]
mu_vs[2] <- data_vs[2]
sd1_vs <- data_vs[3]
sd2_vs <- data_vs[4]
rho_vs <- data_vs[5]
lambda_vs <- data_vs[6]
tot_density_vs <- 0
# calculates the bivariate density for each point
for (i in seq(seq_along(params_vs))) {
Q_vs <- (params_vs[i] - mu_vs[1])^2 / sd1_vs^2 + (fixed_y_coord[i] - mu_vs[2])^2 / sd2_vs^2 -
2 * rho_vs * (params_vs[i] - mu_vs[1]) * (fixed_y_coord[i] - mu_vs[2]) / (sd1_vs * sd2_vs)
density_vs <- 1 / (2 * pi * sd1_vs * sd2_vs * sqrt(1 - rho_vs^2)) * exp(-Q_vs / (2 * (1 -
rho_vs^2)))
tot_density_vs <- tot_density_vs + density_vs
tot_density_vs
}
segment_vs <- add_start_end(toy_trj_vs)
dist_vs <- points_distance(segment_vs)
trj_vs <- sum(dist_vs$distance, na.rm = TRUE)
# total penalty
tot_vs <- lambda_vs * (trj_vs / 10) + (1 - lambda_vs) * (tot_density_vs / 0.1591549)
# return(c(trj_vs, tot_vs))
return(tot_vs)
}
}
trj_vector_penalty_fun <- MAKE_trj_vector_penalty_fun(c(0, 0, 1, 1, 0, 0.5), y_coord)
trj_vector_penalty_fun(x_coord)
#> [1] 0.6677013
initial_vector <- c(0, 0, 0)
optim(par = initial_vector, trj_vector_penalty_fun)$par
结果应该和以前一样，但结果不同：
#> [1] -34.933333 42.885185 9.392593
我检查了这两个函数的所有输出值，它们看起来完全一样。
使用调试器，我在这两种情况下进入 optim()，输出看起来一样，除了一点：
optim(par = c(0, 0, 0), trj_points_penalty_fun)$par 给function = 142
optim(par = initial_vector, trj_vector_penalty_fun)$par 给function = 46
由reprex package (v2.0.0) 于 2021-05-12 创建
【问题讨论】：
-
你能把你的例子减少至少 50% 吗？这可能会帮助您回答自己的问题，但如果没有，其他人可能更愿意看看。
-
我试过了，但我需要两个函数（轨迹函数和二元密度函数）来重现错误。我不知道该剪什么。我可以使段更短，或者不使用构造函数，但这只会节省几行。
标签： r optimization