R中绘制不可导函数时出现异常行为的原因咨询

问题重现 #

先创建数据与分段函数：

v1 <- c(3.323,4.0034,4.6584,5.2826)
v2 <- c(1200,1300,1400,1500)
data <- cbind(v1, v2)

fun <- function(x){
  ans <- NA
  ans[x >= data[1,1] & x <= data[2,1]] <- 100*(x-data[1,1])/(data[2,1]-data[1,1])+data[1,2]
  ans[x >= data[2,1] & x <= data[3,1]] <- 100*(x-data[2,1])/(data[3,1]-data[2,1])+data[2,2]
  ans[x >= data[3,1] & x <= data[4,1]] <- 100*(x-data[3,1])/(data[4,1]-data[3,1])+data[3,2]
  ans
}

单点求值结果符合预期：

> fun(4.3309)
[1] 1350
> fun(4.0034)
[1] 1300

但执行绘图代码时出现警告，且结果不符：

plot(fun, type = 'l', xlim = c(3.323,5.2826))

警告信息：

1: In ans[x >= data[1, 1] & x <= data[2, 1]] <- 100 * (x - data[1, :
number of items to replace is not a multiple of replacement length
2: In ans[x >= data[2, 1] & x <= data[3, 1]] <- 100 * (x - data[2, :
number of items to replace is not a multiple of replacement length
3: In ans[x >= data[3, 1] & x <= data[4, 1]] <- 100 * (x - data[3, :
number of items to replace is not a multiple of replacement length

问题原因 #

plot函数调用自定义函数时，会自动传入一组向量值（而非单个数值）来生成绘图点。原函数中ans初始化为单个NA，当x是向量时，赋值操作的左右两边长度不匹配，导致警告和错误的计算结果。

解决方案 #

修改函数使其支持向量输入，初始化ans为与输入x长度一致的向量，再按区间赋值：

fun <- function(x){
  # 初始化与x同长度的向量
  ans <- numeric(length(x))
  # 第一段区间赋值
  idx1 <- x >= data[1,1] & x <= data[2,1]
  ans[idx1] <- 100*(x[idx1]-data[1,1])/(data[2,1]-data[1,1])+data[1,2]
  # 第二段区间赋值
  idx2 <- x >= data[2,1] & x <= data[3,1]
  ans[idx2] <- 100*(x[idx2]-data[2,1])/(data[3,1]-data[2,1])+data[2,2]
  # 第三段区间赋值
  idx3 <- x >= data[3,1] & x <= data[4,1]
  ans[idx3] <- 100*(x[idx3]-data[3,1])/(data[4,1]-data[3,1])+data[3,2]
  # 超出区间的设为NA
  ans[!(idx1 | idx2 | idx3)] <- NA
  ans
}

也可以用ifelse实现更简洁的向量化逻辑：

fun <- function(x){
  ifelse(x >= data[1,1] & x <= data[2,1],
         100*(x-data[1,1])/(data[2,1]-data[1,1])+data[1,2],
         ifelse(x >= data[2,1] & x <= data[3,1],
                100*(x-data[2,1])/(data[3,1]-data[2,1])+data[2,2],
                ifelse(x >= data[3,1] & x <= data[4,1],
                       100*(x-data[3,1])/(data[4,1]-data[3,1])+data[3,2],
                       NA)))
}

修改后再执行plot(fun, type='l', xlim=c(3.323,5.2826))，无警告且绘图结果与单点求值一致。

内容的提问来源于stack exchange，提问作者Edward Erisman