同じコードを持つ2つのRスクリプト：1つは結果、もう1つはスローエラーです

Question

2つの同一のRスクリプトがあります。しかし、私がそれぞれ実行すると、1つが動作しますが、もう1つはエラーをスローします。私は、それぞれ異なるRセッションで実行し、それらを実行する前にメモリをクリアします。

なぜこれが当てはまるのかについての回答は見つかりませんでした。

おそらく、他の誰かがこの問題を抱えていて、なぜそれが発生しているのか知っていますか？

スクリプトの作業

library(investr) 
library(mgcv) 
library(rootSolve) 
library(scam) 

inv.predict2 <- function(object, y, x.name, interval = FALSE, 
        lower, upper, level = 0.95,...) { 
.fun1 <- function(x) { 
    predFit(object, newdata = setNames(data.frame(x), x.name)) - y 
} 
.fun2 <- function(x) { 
    predFit(object, newdata = setNames(data.frame(x), x.name), 
    interval = "confidence")[, "upr"] - y 
    } 
.fun3 <- function(x) { 
    predFit(object, newdata = setNames(data.frame(x), x.name), 
     interval = "confidence")[, "lwr"] - y 
} 
x0.est <- uniroot.all(lower = lower, upper = upper, ..., f = .fun1) 
     res <- if (interval) { 
      lwr <- uniroot.all(lower = lower, upper = x0.est, ..., f = .fun2) 
      upr <- uniroot.all(lower = x0.est, upper = upper, ..., f = .fun3) 
      lwr <- min(c(lwr, upr)) 
      upr <- max(c(lwr, upr)) 
c("estimate" = x0.est, "lower" = lwr, "upper" = upr) 
} else { 
x0.est 
} 
res 
} 

predFit.gam <- function(object, newdata, type = c("link", "response"), 
        interval = c("none", "confidence", "prediction"), 
        level = 0.95, ...) { 
type <- match.arg(type) 
interval <- match.arg(interval) 
res <- if (interval == "none") { 
predict.gam(object, newdata = newdata, type = type, ...) 
} else if (interval == "confidence") { 
pred <- predict.gam(object, newdata = newdata, se.fit = TRUE, type = "link", 
       ...) 
out <- cbind("fit" = pred$fit, 
      "lwr" = pred$fit - pred$se.fit * stats::qnorm((level + 1)/2), 
      "upr" = pred$fit + pred$se.fit * stats::qnorm((level + 1)/2)) 
if (type == "response") { 
    out <- apply(out, MARGIN = 2, FUN = function(x) { 
    stats::family(object)$linkinv(x) 
    }) 
} 
out 
} else { 
stop("Prediction intervals are currently not supported for GAMs.") 
} 
res 
} 

predFit.scam <- function(object, newdata, type = c("link", "response"), 
        interval = c("none", "confidence", "prediction"), 
        level = 0.95, ...) { 
    type <- match.arg(type) 
    interval <- match.arg(interval) 
    res <- if (interval == "none") { 
predict.scam(object, newdata = newdata, type = type, ...) 
} else if (interval == "confidence") { 
pred <- predict.scam(object, newdata = newdata, se.fit = TRUE, type = "link", 
       ...) 
out <- cbind("fit" = pred$fit, 
      "lwr" = pred$fit - pred$se.fit * stats::qnorm((level + 1)/2), 
      "upr" = pred$fit + pred$se.fit * stats::qnorm((level + 1)/2)) 
if (type == "response") { 
    out <- apply(out, MARGIN = 2, FUN = function(x) { 
    stats::family(object)$linkinv(x) 
    }) 
} 
out 
} else { 
stop("Prediction intervals are currently not supported for SCAMs.") 
} 
res 
} 


ptm <- proc.time() 

set.seed(1) 

# Rprof() 

K <- 1 
N <- 100 
Hstar <- 10 

perms <- 10000 

specs <- 1:N 

pop <- array(dim = c(c(perms, N), K)) 

haps <- as.character(1:Hstar) 

probs <- rep(1/Hstar, Hstar) 

x <- c(1:3) 
y <- c(3:5) 

for(j in 1:perms){ 
    for(i in 1:K){ 
     if(i == 1){ 
     pop[j, specs, i] <- sample(haps, size = N, replace = TRUE, prob = probs) 
    } 
     else{ 
      pop[j ,, 1] <- sample(haps[x], size = N, replace = TRUE, prob = probs[x]) 
      pop[j ,, 2] <- sample(haps[y], size = N, replace = TRUE, prob = probs[y]) 
     } 
} 
} 

HAC.mat <- array(dim = c(c(perms, N), K)) 

for(k in specs){ 
    for(j in 1:perms){ 
     for(i in 1:K){ 
      ind.index <- sample(specs, size = k, replace = FALSE) 
      hap.plot <- pop[sample(1:nrow(pop), size = 1, replace = TRUE), ind.index, sample(1:K, size = 1, replace = TRUE)] 
      HAC.mat[j, k, i] <- length(unique(hap.plot)) 
    } 
} 
} 

means <- apply(HAC.mat, MARGIN = 2, mean) 
lower <- apply(HAC.mat, MARGIN = 2, function(x) quantile(x, 0.025)) 
upper <- apply(HAC.mat, MARGIN = 2, function(x) quantile(x, 0.975)) 

par(mfrow = c(1, 2)) 

plot(specs, means, type = "n", xlab = "Specimens sampled", ylab = "Unique haplotypes", ylim = c(1, Hstar)) 
polygon(x = c(specs, rev(specs)), y = c(lower, rev(upper)), col = "gray") 
lines(specs, means, lwd = 2) 
HAC.bar <- barplot(N*probs, xlab = "Unique haplotypes", ylab = "Specimens sampled", names.arg = 1:Hstar) 

# summaryRprof() 

proc.time() - ptm 

d <- data.frame(specs, means) 

HAC.tp <- gam(means ~ s(specs, bs = "tp", k = 20), optimizer = c("outer", "bfgs"), data = d) # thin plate spline 

HAC.tp <- inv.predict2(HAC.tp, y = Hstar, x.name = "specs", interval = TRUE, lower = 1, upper = 1000000) 
HAC.tp 

非作業スクリプト

library(investr) 
library(mgcv) 
library(rootSolve) 
library(scam) 

inv.predict2 <- function(object, y, x.name, interval = FALSE, 
        lower, upper, level = 0.95,...) { 
.fun1 <- function(x) { 
predFit(object, newdata = setNames(data.frame(x), x.name)) - y 
} 
.fun2 <- function(x) { 
predFit(object, newdata = setNames(data.frame(x), x.name), 
     interval = "confidence")[, "upr"] - y 
} 
.fun3 <- function(x) { 
predFit(object, newdata = setNames(data.frame(x), x.name), 
     interval = "confidence")[, "lwr"] - y 
} 
x0.est <- uniroot.all(lower = lower, upper = upper, ..., f = .fun1) 
res <- if (interval) { 
lwr <- uniroot.all(lower = lower, upper = x0.est, ..., f = .fun2) 
upr <- uniroot(lower = x0.est, upper = upper, ..., f = .fun3) 
lwr <- min(c(lwr, upr)) 
upr <- max(c(lwr, upr)) 
c("estimate" = x0.est, "lower" = lwr, "upper" = upr) 
} else { 
x0.est 
} 
res 
} 

predFit.gam <- function(object, newdata, type = c("link", "response"), 
        interval = c("none", "confidence", "prediction"), 
        level = 0.95, ...) { 
type <- match.arg(type) 
interval <- match.arg(interval) 
res <- if (interval == "none") { 
predict.gam(object, newdata = newdata, type = type, ...) 
} else if (interval == "confidence") { 
pred <- predict.gam(object, newdata = newdata, se.fit = TRUE, type = "link", 
       ...) 
out <- cbind("fit" = pred$fit, 
      "lwr" = pred$fit - pred$se.fit * stats::qnorm((level + 1)/2), 
      "upr" = pred$fit + pred$se.fit * stats::qnorm((level + 1)/2)) 
if (type == "response") { 
    out <- apply(out, MARGIN = 2, FUN = function(x) { 
    stats::family(object)$linkinv(x) 
    }) 
} 
out 
} else { 
stop("Prediction intervals are currently not supported for GAMs.") 
} 
res 
} 

predFit.scam <- function(object, newdata, type = c("link", "response"), 
        interval = c("none", "confidence", "prediction"), 
        level = 0.95, ...) { 
    type <- match.arg(type) 
    interval <- match.arg(interval) 
    res <- if (interval == "none") { 
predict.scam(object, newdata = newdata, type = type, ...) 
} else if (interval == "confidence") { 
pred <- predict.scam(object, newdata = newdata, se.fit = TRUE, type = "link", 
       ...) 
out <- cbind("fit" = pred$fit, 
      "lwr" = pred$fit - pred$se.fit * stats::qnorm((level + 1)/2), 
      "upr" = pred$fit + pred$se.fit * stats::qnorm((level + 1)/2)) 
if (type == "response") { 
    out <- apply(out, MARGIN = 2, FUN = function(x) { 
    stats::family(object)$linkinv(x) 
    }) 
} 
out 
} else { 
stop("Prediction intervals are currently not supported for SCAMs.") 
} 
res 
} 

ptm <- proc.time() 

set.seed(1) 

# Rprof() 

K <- 1 

N <- 100 

Hstar <- 10 

perms <- 10000 

specs <- 1:N 

pop <- array(dim = c(c(perms, N), K)) 

haps <- as.character(1:Hstar) 

probs <- rep(1/Hstar, Hstar) 

s1 <- c(1:6) 

s2 <- c(7:10) 

for(j in 1:perms){ 
    for(i in 1:K){ 
     if(i == 1){ 
      pop[j, specs, i] <- sample(haps, size = N, replace = TRUE, prob = probs) 
    } 
     else{ 
      pop[j ,, 1] <- sample(haps[s1], size = N, replace = TRUE, prob = probs[s1]) 
      pop[j ,, 2] <- sample(haps[s2], size = N, replace = TRUE, prob = probs[s2]) 
     } 
} 
} 

HAC.mat <- array(dim = c(c(perms, N), K)) 

for(k in specs){ 
    for(j in 1:perms){ 
     for(i in 1:K){ 
      ind.index <- sample(specs, size = k, replace = FALSE) 
      hap.plot <- pop[sample(1:nrow(pop), size = 1, replace = TRUE), ind.index, sample(1:K, size = 1, replace = TRUE)] 
      HAC.mat[j, k, i] <- length(unique(hap.plot)) 
    } 
} 
} 

means <- apply(HAC.mat, MARGIN = 2, mean) 
lower <- apply(HAC.mat, MARGIN = 2, function(x) quantile(x, 0.025)) 
upper <- apply(HAC.mat, MARGIN = 2, function(x) quantile(x, 0.975)) 

par(mfrow = c(1, 2)) 

plot(specs, means, type = "n", xlab = "Specimens sampled", ylab = "Unique haplotypes", ylim = c(1, Hstar)) 
polygon(x = c(specs, rev(specs)), y = c(lower, rev(upper)), col = "gray") 
lines(specs, means, lwd = 2) 
HAC.bar <- barplot(N*probs, xlab = "Unique haplotypes", ylab = "Specimens sampled", names.arg = 1:Hstar) 

# summaryRprof() 

proc.time() - ptm 

d <- data.frame(specs, means) 

HAC.tp <- gam(means ~ s(specs, bs = "tp", k = 20), optimizer = c("outer", "bfgs"), data = d) # thin plate spline 
summary(HAC.tp) 
plot(HAC.tp) 

HAC.tp <- inv.predict2(HAC.tp, y = Hstar, x.name = "specs", interval = TRUE, lower = 1, upper = 1000000) 
HAC.tp 

Answer 1

これらの2つのスクリプトは同じではありません。

22行目の作業スクリプトでは、uniroot.all()があり、動作していないスクリプトにはuniroot()があります。私は非作業をuniroot.all()に変更し、エラーなしで完了しました（しかし、結果が正しいかどうかは全く考えられません）。

も ​​- その有用な情報ならば、ライン102の周りに、あなたは（作業）の違いがあります。

x <- c(1:3) 
y <- c(3:5) 

（非稼働）を：

s1 <- c(1:6) 

s2 <- c(7:10) 

この違いは、おそらくそれほど重要ではありません私はあなたに知らせると思っただけです。

最後に、今後、BeyondCompare、またはAtomの使用をお勧めします。 （これらの2つのプログラムは、テキスト（またはファイルやフォルダ）を比較して簡単に見つけることができます違い。

Answer 2

たぶんデバッグ機能を探索するためにこれらのいずれかを試してみてください？

traceback() 
    browser() 

http://adv-r.had.co.nz/Exceptions-Debugging.html