将R中的两个相似函数矢量化为一个

Today, I came across a problem: two almost identical functions work as intended before vectorisation, but after it, one works fine, and another one returns an error.

今天,我遇到了一个问题:两个几乎相同的函数在向量化之前按预期工作,但在它之后,一个工作正常,另一个返回错误。

I am examining the robustness of various estimators with respect to different transformations of residuals and aggregating functions. Quantile Regression and Least Median of Squares are particular cases of what I am doing.

我正在研究各种估计器在残差和聚合函数的不同变换方面的鲁棒性。分位数回归和最小正方形是我正在做的事情的特例。

So I wrote the following code to see how the Least Trimean of Squares is going to work and found out that it works fine if model parameters are supplied as different arguments, but fails if they come in a vector. For instance, I need the first function for plotting (it is convenient to use outer(...) to get a value matrix for persp or just supply f(x, y) to persp3d from library(rgl), but the second one for estimation (R optimisers are expecting a vector of inputs as the first argument over which the minimisation is going to be done).

因此,我编写了以下代码,以了解Squares的Least Trimean是如何工作的,并发现如果模型参数作为不同的参数提供,它可以正常工作,但如果它们出现在向量中则会失败。例如,我需要第一个绘图功能(使用outer(...)来获取persp的值矩阵或者只是从库(rgl)向persp3d提供f(x,y)是方便的,但是第二个用于估计(R优化器期望输入向量作为最小化将要进行的第一个参数)。

MWE:

set.seed(105)
N <-  204
x <- rlnorm(N)
y <- 1 + x + rnorm(N)*sqrt(.1+.2*x+.3*x^2)
# A simple linear model with heteroskedastic errors

resfun <- function(x) return(x^2)
# Going to minimise a function of squared residuals...
distfun <- function(x) return(mean(quantile(x, c(0.25, 0.5, 0.5, 0.75)))) 
# ...which in this case is the trimean

penalty <- function(theta0, theta1) {
  r <- y - theta0 - theta1*x
  return(distfun(resfun(r)))
}

pen2 <- function(theta) {
  r <- y - theta[1] - theta[2]*x
  return(distfun(resfun(r)))
}

penalty(1, 1) # 0.5352602
pen2(c(1, 1)) # 0.5352602

vpenalty <- Vectorize(penalty)
vpen2 <- Vectorize(pen2)

vpenalty(1, 1) # 0.5352602
vpen2(c(1, 1))

Error in quantile.default(x, c(0.25, 0.5, 0.5, 0.75)) : 
missing values and NaN's not allowed if 'na.rm' is FALSE

Why does vpen2, being vectorised pen2, choke even on a single input?

为什么vpen2,即矢量化笔2,甚至在单个输入上窒息?

1 个解决方案

#1

As jogo pointed out, vpen2 reads the elements of the input vector and tries to take the first one. The right way to go is to use something like

正如jogo所指出的,vpen2读取输入向量的元素并尝试获取第一个元素。正确的方法是使用类似的东西

a <- matrix(..., ncol=2)
apply(a, 1, pen2)

This will return a vector of values from vpar2 evaluated at each row of the matrix.

这将返回在矩阵的每一行评估的vpar2的值向量。

#1