1. 首頁
  2. 問答
  3. 具有相同代碼的兩個R腳本：一個給出結果，另一個引發錯誤

具有相同代碼的兩個R腳本：一個給出結果，另一個引發錯誤

2017-08-25 50 views
0

我有兩個相同的R腳本;然而，當我運行每一個，一個工作，但另一個拋出一個錯誤。我運行它們之前，在不同的R會話中運行它們並清除內存。具有相同代碼的兩個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

來源

2017-08-25 Jarrett Phillips

+1

將腳本粘貼到此處很有幫助。 – Roshan

+0

只是猜測而不看代碼。也許是隨機抽樣/選擇引發的錯誤？嘗試set.seed（）。 –

+0

查看腳本可能有助於提供解決方案。 – Sagar

回答

0

這兩個腳本是不相同的。

在第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。 （我想我更喜歡Atom，但它已經有一段時間了，因爲我們必須在這裏使用BC）。這兩個程序（或許多其他類似的替代品）可以幫助您比較文本（或文件或文件夾）以輕鬆找到差異。

來源

2017-08-25 20:11:26 lukehawk

+0

謝謝 - 不知道這是怎麼發生的。完全錯過了......我知道第102行的區別，但它不會影響結果。我以前從未聽說beyondCompare和Atom ......我一定會檢查它們！ –

+0

沒有問題！很高興我終於可以回答一個問題，而不是總是問他們！ – lukehawk

+0

我總是問問題！ –

相關問題

 相關問題