如何將參數傳遞到R中的mapply函數？

Question

背景。 我正在閱讀the paper並試圖找到(tau1*, tau2*) = arg max P_D(tau1, tau2)（等式（30））。在論文（第6頁，表1）中，您可以看到作者所獲得的結果（專欄 - 主席 - 瓦斯尼規則）。我用手在[1，15]範圍內改變了初始參數tau1，tau2，我的結果接近原始結果。

該圖顯示了初始參數爲tau1=tau2=1（藍線）和tau1=tau2=15（紅線）並與「Chair-Varshney規則」（黑點）比較時的結果。

我的代碼如下。

fun_PD <- function(par, alpha, N){ 

t1 <- par[[1]]; t2 <- par[[2]] 
lambdab <- 10 
lambdac <- c(0.625, 0.625) 
sigma2_w <- 10 
p<-c(); q<-c() 

# Compute P-values, complementary CDF 
p[1]<- 1 - pnorm((t1 - lambdab - lambdac[1])/sqrt(sigma2_w + lambdab + lambdac[1])) # (5) 
p[2]<- 1 - pnorm((t2 - lambdab - lambdac[2])/sqrt(sigma2_w + lambdab + lambdac[2])) # (6) 

q[1] <- 1 - pnorm((t1 - lambdab)/sqrt(sigma2_w + lambdab)) # (7) 
q[2] <- 1 - pnorm((t2 - lambdab)/sqrt(sigma2_w + lambdab)) # (8) 

Q00 <- (1-q[1])*(1-q[2]); Q01 <- (1-q[1])*q[2] # page 4 
Q10 <- q[1]*(1-q[2]);  Q11 <- q[1]*q[2] 

P00 <- (1-p[1])*(1-p[2]); P01 <- (1-p[1])*p[2] # page 5 
P10 <- p[1]*(1-p[2]);  P11 <- p[1]*p[2] 

C <- c(log((P10*Q00)/(P00*Q10)), log((P01*Q00)/(P00*Q01))) # (13) 

mu0 <- N * (C[1]*q[1] + C[2]*q[2]) # (14) 
mu1 <- N * (C[1]*p[1] + C[2]*p[2]) # (16) 

sigma2_0 <- N * (C[1]^2*q[1]*(1-q[1]) + C[2]^2*q[2]*(1-q[2])) # (15) 
sigma2_1 <- N * (C[1]^2*q[1]*(1-q[1]) + C[2]^2*q[2]*(1-q[2])) # (17) 

sigma0 <- sqrt(sigma2_0) 
sigma1 <- sqrt(sigma2_1) 

#Compute critical values, inverse of the CCDF 
PA <- qnorm(alpha, lower.tail=FALSE) 

gamma <- sigma0 * PA + mu0     # (20) 
out <- 1 - pnorm((gamma - mu1)/sigma1) # (30) 

return(out) 
} # fun_PD 
########################################################################### 

dfb <- data.frame(a=c(0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5), 
        r=c(.249, .4898, .6273, .7738, .8556, .9076, .9424)) 

df <- data.frame() 
a <- seq(0,1,0.05) 
n <- length(a) 
for(i in 1:n) { 
tau_optimal <- optim(par=c(t1=1,t2=1),   # parameter 
        fn=fun_PD, 
        control=list(fnscale=-1), # maximization 
        method="CG", 
        alpha = a[i],    # const 
        N = 100)     # const 
df = rbind(df, c(tau_optimal$par[1], tau_optimal$par[2], a[i], tau_optimal$value)) 
} 
colnames(df) <- c("tau1", "tau2", "alpha", "P_d") 
df 

一些模擬後，我understud該函數fun_P_D可以有一些地方的最低高度和極大值，我曾嘗試使用圖形approuch從R-User-guide檢測功能的本地最低高度和極大值：

該 馬塞洛的更新答案

編輯2.後：

fun_PDtest <- function(x, y){ 
mapply(fun_PD, x, y, MoreArgs = list(N=100, alpha=0.1)) 
} 
x<-(1:10); y<-c(1:10) 
fun_PDtest(x,y) 
# Error in (function (par, alpha, N) : unused argument (dots[[2]][[1]]) 

我的問題是：如何將向量x,y轉換成mapply函數？

Answer 1

outer擴展了2個向量，並期望函數帶有2個相同大小的向量。而不是重寫fun_PD採取向量，您可以使用mapply並在fun_PDtest中調用原始函數。您還可以創建接收optmin

中使用的矢量完整代碼的函數：

#Rewrite function to use x, y instead of receiving a vector 
fun_PD <- function(x , y, alpha, N) { 

    t1<-y 
    t2<-x 

    N<-100 
    alpha<-0.1 
    lambdab <- 10 
    lambdac <- c(0.625, 0.625) 
    sigma2_w <- 10 
    p<-c(); q<-c() 

    # Compute P-values, complementary CDF 
    p[1]<- 1 - pnorm((t1 - lambdab - lambdac[1])/sqrt(sigma2_w + lambdab + lambdac[1])) # (5) 
    p[2]<- 1 - pnorm((t2 - lambdab - lambdac[2])/sqrt(sigma2_w + lambdab + lambdac[2])) # (6) 

    q[1] <- 1 - pnorm((t1 - lambdab)/sqrt(sigma2_w + lambdab)) # (7) 
    q[2] <- 1 - pnorm((t2 - lambdab)/sqrt(sigma2_w + lambdab)) # (8) 

    Q00 <- (1-q[1])*(1-q[2]); Q01 <- (1-q[1])*q[2] # page 4 
    Q10 <- q[1]*(1-q[2]);  Q11 <- q[1]*q[2] 

    P00 <- (1-p[1])*(1-p[2]); P01 <- (1-p[1])*p[2] # page 5 
    P10 <- p[1]*(1-p[2]);  P11 <- p[1]*p[2] 

    C <- c(log((P10*Q00)/(P00*Q10)), log((P01*Q00)/(P00*Q01))) # (13) 

    mu0 <- N * (C[1]*q[1] + C[2]*q[2]) # (14) 
    mu1 <- N * (C[1]*p[1] + C[2]*p[2]) # (16) 

    sigma2_0 <- N * (C[1]^2*q[1]*(1-q[1]) + C[2]^2*q[2]*(1-q[2])) # (15) 
    sigma2_1 <- N * (C[1]^2*q[1]*(1-q[1]) + C[2]^2*q[2]*(1-q[2])) # (17) 

    sigma0 <- sqrt(sigma2_0) 
    sigma1 <- sqrt(sigma2_1) 

    #Compute critical values, inverse of the CCDF 
    PA <- qnorm(alpha, lower.tail=FALSE) 

    gamma <- sigma0 * PA + mu0     # (20) 
    out <- 1 - pnorm((gamma - mu1)/sigma1) # (30) 

    return(out) 

} 

x<-seq(1,15, len=50) 
y<-seq(1,15, len=50) 

# then I rewrite my function without passing alpha and N 

fun_PDimage <- function(x, y){ 

mapply(fun_PD,x,y, MoreArgs = list(N=100, alpha=0.1)) 
    # the body is the same as in fun_PD(par, alpha, N) 
} # fun_PDimage 

z <-outer(x, y, fun_PDimage) # errors are here 

# Rewrite function for use in optim 
fun_PDoptim <- function(v){ 

    x<-v[1] 
    y<-v[2] 

    fun_PD(x, y, 0.1, 100) 
} # fun_PDoptim 

#Create the image 
image(x,y,z, col=heat.colors(100)) 
contour(x,y,z,add=T) 

# Find the max using optmin 
res<-optim(c(2,2),fun_PDoptim, control = list(fnscale=-1)) 
print(res$par) 

#Add Point to image 
points(res$par[1], res$par[2],pch=3) 

下面是結果： 點，其中的功能有最大：

> print(res$par) 
[1] 12.20753 12.20559 

圖片：