「在...％*％...中出錯：非適形參數」在迴歸中使用自己的函數

Question

我有一個函數Q（x | delta）：R^n - > R我想擬合非線性分位數迴歸。函數Q（。）使用了一些矩陣運算，不使用它會非常複雜。問題是，當在公式參數中使用的函數中存在矩陣運算時，似乎nlrq（非線性分位數迴歸）和nls（非線性迴歸）不起作用。爲了說明，考慮在我不使用矩陣運算但不使用矩陣運算的情況下我可以在nlrq和nls函數的公式變量中使用的更簡單的函數F（x1，x2 | a，b，c）在用矩陣運算寫入時在公式參數中工作。

library('quantreg') 

    ## Generating the data 
    x1<- rnorm(200) 
    x2<- rnorm(200) 
    y<- 1+3*sin(x1)+2*cos(x2) +rnorm(200) 
    Dat<- data.frame(y,x1,x2) 

    ## The function F1 without matrix operation 
    F1<- function(x_1, x_2, a, b,c){a+b*sin(x_1)+c*cos(x_2)} 

    ## The function F2 with matrix operation 
    F2<- function(x_1, x_2, a, b,c){t(c(1,sin(x_1),cos(x_2)))%*%c(a,b,c)} 

    ## Both functions work perfectly 
    F1(x_1=3, x_2=2, a=1, b=3,c=2) 
    F2(x_1=3, x_2=2, a=1, b=3,c=2) 

    ## But only F1 can be estimated by nls and nlrq 
    nls_1<-nls(y ~ F1(x_1 = x1, x_2 = x2, a = 1, b, c), 
       data = Dat, start = list(b = 3, c = 2)) 

    nlrq_1<-nlrq(y ~ F1(x_1 = x1, x_2 = x2, a = 1, b, c), 
       data = Dat, start = list(b = 3, c = 2), tau = 0.9) 

    ## When F2 is used in the formula argument an error happens 
    nls_2<-nls(y ~ F2(x_1 = x1, x_2 = x2, a = 1, b, c), 
       data = Dat, start = list(b = 3, c = 2)) 

    nlrq_2<-nlrq(y ~ F2(x_1 = x1, x_2 = x2, a = 1, b, c), 
       data = Dat, start = list(b = 3, c = 2), tau = 0.9) 

的錯誤是Error in t(c(1, sin(x_1), cos(x_2))) %*% c(a, b, c) : non-conformable arguments。我相信如果有人設法通過使用nls和nlrq的矩陣運算來估計F2，我可以在我的其他函數中使用相同的解決方案。

Dat的大小是200x3。

非常感謝。

Answer 1

您的功能F2()不適用於向量參數x_1,x_2，...因爲c(...)只構造一個長向量（不是矩陣）。
參見：

F1(x_1=c(3,5), x_2=c(2,4), a=1, b=3,c=2) 
F2(x_1=c(3,5), x_2=c(2,4), a=1, b=3,c=2) 

結果：

#> F1(x_1=c(3,5), x_2=c(2,4), a=1, b=3,c=2) 
#[1] 0.5910664 -3.1840601 
#> F2(x_1=c(3,5), x_2=c(2,4), a=1, b=3,c=2) 
#error in t(c(1, sin(x_1), cos(x_2))) %*% c(a, b, c) : ... 

功能nls()和nlrq()被（從數據幀Dat即列）發送向量的函數F2()（相應F1()）。

下面是一些F2()矢量定義：

# other definitions for F2() 
F2 <- function(x_1, x_2, a, b,c) cbind(1,sin(x_1),cos(x_2)) %*% c(a,b,c) 
F2(x_1=c(3,5), x_2=c(2,4), a=1, b=3,c=2) 

F2 <- function(x_1, x_2, a, b,c) t(rbind(1,sin(x_1),cos(x_2))) %*% c(a,b,c) 
F2(x_1=c(3,5), x_2=c(2,4), a=1, b=3,c=2) 

F2 <- function(x_1, x_2, a, b,c) colSums(rbind(1,sin(x_1),cos(x_2)) * c(a,b,c)) 
F2(x_1=c(3,5), x_2=c(2,4), a=1, b=3,c=2) 

F2 <- function(x_1, x_2, a, b,c) crossprod(rbind(1,sin(x_1),cos(x_2)), c(a,b,c)) 
F2(x_1=c(3,5), x_2=c(2,4), a=1, b=3,c=2) 

F2 <- function(x_1, x_2, a, b,c) tcrossprod(c(a,b,c), cbind(1,sin(x_1),cos(x_2))) 
F2(x_1=c(3,5), x_2=c(2,4), a=1, b=3,c=2) 

Answer 2

可以回落到這個通用的優化功能。 R中通常的默認值是optim，但還有很多其他值。

下面是最小二乘迴歸的情況。損失函數是平方殘差的總和。我已經重寫了你的F2函數，以便它適用於向量參數。

sumsq <- function(beta) 
{ 
    F2 <- function(x1, x2, a, b, c) 
    { 
     cbind(1, sin(x1), cos(x2)) %*% c(a, b, c) 
    } 
    yhat <- F2(Dat$x1, Dat$x2, beta[1], beta[2], beta[3]) 
    sum((Dat$y - yhat)^2) 
} 

beta0 <- c(mean(Dat$y), 1, 1) 

optim(beta0, sumsq, method="BFGS") 

#initial value 731.387431 
#final value 220.265745 
#converged 
#$par 
#[1] 0.8879371 3.0211286 2.1639280 
# 
#$value 
#[1] 220.2657 
# 
#$counts 
#function gradient 
#  25  7 
# 
#$convergence 
#[1] 0 
# 
#$message 
#NULL 

這裏，optim返回一個列表與多個組件。組件par是最小化組件value中的方形殘差總和的迴歸係數的值。

如果與nls的結果進行比較，您可以看到估計的係數大致相等。

nls(y ~ F1(x_1=x1, x_2=x2, a=1, b, c), 
      data=Dat, start=list(b=3, c=2)) 

Nonlinear regression model 
    model: y ~ F1(x_1 = x1, x_2 = x2, a = 1, b, c) 
    data: Dat 
    b  c 
3.026 2.041 
residual sum-of-squares: 221 

Number of iterations to convergence: 1 
Achieved convergence tolerance: 7.823e-10 

你可以做一些類似的分位數迴歸，但那會更復雜。

Answer 3

基於其他答案，我已經發現問題是使用c()函數在F2中構建向量。當我使用rbind()代替時，估計與nls()和nlrq()完美結合。

接下來我將顯示F2的更正版本。

## Changing c() for rbind() 
    F2<- function(x_1, x_2, a, b,c){t(rbind(1,sin(x_1),cos(x_2)))%*%rbind(a,b,c)} 

    ## Now nls() and nlrq() work properly 
    nls_2<-nls(y ~ F2(x_1 = x1, x_2 = x2, a = 1, b, c), 
     data = Dat, start = list(b = 3, c = 2)) 

    nlrq_2<-nlrq(y ~ F2(x_1 = x1, x_2 = x2, a = 1, b, c), 
     data = Dat, start = list(b = 3, c = 2), tau = 0.9) 

注意，在nls_2和估計nlrq_2與來自nls_1和nlrq_1的那些一致。

非常感謝您的幫助。