2012-08-02 50 views
4

這是我遇到的一個非常奇怪的情況。基本上,我試圖將累積分佈函數擬合到我的數據的G函數中。完成後,我想繪製模型和原始數據,並將其輸出爲PDF。我會讓代碼解釋(簡單的複製和粘貼):導致對象消失的函數

library(spatstat) 

data(swedishpines) 

mydata <- swedishpines 

mydata.Gest <- Gest(mydata) 

Gvalues <- mydata.Gest$rs 

count <- (which(Gvalues == 1))[1] 

new_r <- seq(1/count, length(Gvalues)/count, by = 1/count) 

GvsR_dataframe <- data.frame(G <- Gvalues, R <- new_r) 

themodel <- suppressWarnings(nls(G ~ pnorm(R, mean, sd), data = GvsR_dataframe, start = list(mean=0.4, sd=0.2), trace = FALSE)) 

pdf(file = "ModelPlot.pdf") 

plot(mydata.Gest, cbind(rs, theo) ~ new_r, lty = c(1, 2), col = c("black", "red"), xlim = c(0, max(new_r)), ylim = c(0,1), main = paste("Model-fitting for G Function \n Mean = ",as.numeric(coef(themodel)[1]),"\n Standard Deviation = ",as.numeric(coef(themodel)[2]), sep=''), ylab = "G(r)", xlab = "Distance Between Particles (r)", legend = NULL) 
lines(new_r, predict(themodel), lty = 2, col = "blue") 
legend("bottomright", c("CSR", "Swedish Pines", "Normal Probability \n Density Function"), lty = c(2, 4, 1, 2), col = c("red", "black", "blue"), bg = 'grey', border = 'black') 

graphics.off() 

上面的代碼完美地工作。

現在是奇怪的部分。

當我封裝mydata <- swedishpines作爲函數之後的所有命令,並且導致mydata成爲此函數的輸入時,它不再起作用。以下代碼應該執行代碼的最後一部分,但它不。

library(spatstat) 

data(swedishpines) 

mydata <- swedishpines 

ModelFit <- function(mydata) { 

mydata.Gest <- Gest(mydata) 

Gvalues <- mydata.Gest$rs 

count <- (which(Gvalues == 1))[1] 

new_r <- seq(1/count, length(Gvalues)/count, by = 1/count) 

GvsR_dataframe <- data.frame(G <- Gvalues, R <- new_r) 

themodel <- suppressWarnings(nls(G ~ pnorm(R, mean, sd), data = GvsR_dataframe, start = list(mean=0.4, sd=0.2), trace = FALSE)) 

pdf(file = "ModelPlot.pdf") 

plot(mydata.Gest, cbind(rs, theo) ~ new_r, lty = c(1, 2), col = c("black", "red"), xlim = c(0, max(new_r)), ylim = c(0,1), main = paste("Model-fitting for G Function \n Mean = ",as.numeric(coef(themodel)[1]),"\n Standard Deviation = ",as.numeric(coef(themodel)[2]), sep=''), ylab = "G(r)", xlab = "Distance Between Particles (r)", legend = NULL) 
lines(new_r, predict(themodel), lty = 2, col = "blue") 
legend("bottomright", c("CSR", "Swedish Pines", "Normal Probability \n Density Function"), lty = c(2, 4, 1, 2), col = c("red", "black", "blue"), bg = 'grey', border = 'black') 

graphics.off() 

} 

ModelFit(mydata) 

出現以下錯誤:

Error in eval(expr, envir, enclos) : object 'new_r' not found 

我很困惑。我一直在研究這個問題很長一段時間,只是不能想出解決這個問題的辦法。 PDF輸出,但它是腐敗的,不會打開。我不知道爲什麼new_r'消失',但這樣做會導致所有繪圖操作停止。顯然,new_r是函數ModelFit的本地函數,但它幾乎看起來好像在函數的某些區域也是本地的。

任何幫助將不勝感激。

+0

我無法重現。你確定你的第二個代碼塊中的代碼與你在R中執行的代碼完全相同嗎? – 2012-08-02 19:34:00

+0

我相信這是。我在第二個模塊中改變的是一個函數定義(現在包含第一個模塊的大部分代碼),隨後調用該函數。 – MikeZ 2012-08-02 19:38:56

+0

有趣。我非常感謝任何幫助,因爲我試圖在Linux服務器上以批處理模式運行一個大型程序,而這一小段代碼讓我感到非常悲傷。 – MikeZ 2012-08-02 19:40:10

回答

6

你在那裏做了很多隱含的東西......我會建議更明確地寫東西。

具體而言,mydata.Gest$r <- new_r然後在您的繪圖公式plot(..., cbind(rs, theo) ~ r, ...)中用r代替new_r。這對我行得通。不知道爲什麼它在一個函數之外,而不是在裏面,但依靠plot來看看mydata.Gest的本地範圍以外的new_r是有風險的。

此外,使用=分配的事情列一個數據幀,而不是<-

從一個乾淨的會話:

data.frame(x<-1:10, y<- 1:10) 
ls() 

data.frame(x=1:10, y=1:10) 
ls() 
+0

非常好的解決方案,賈斯汀,謝謝! 但是,我仍然因爲我們沒有針對發生錯誤的事實而拖延。我使用了之前發佈的流程,而且我從來沒有遇到任何問題。我以前的用法和我發佈的用法之間的區別在於,其他用戶都不涉及'nls'。我看到的唯一可能的錯誤是,在使用'new_r'創建'GvsR_dataframe'後,'new_r'對象被'nls'使用後可能會被刪除。我很好奇,想深入瞭解這個問題。 – MikeZ 2012-08-03 13:18:14

+1

我認爲它與'plot.Gest'的寫法有很大關係。它在那個你遇到問題的呼叫中。你可以在劇情調用的其他部分使用'new_r',但公式部分看起來並不超出數據參數的局部範圍。 – Justin 2012-08-03 13:55:06

+0

在使用獨立變量作爲繪圖函數之外的對象的上下文中,我曾經使用過'plot.Gest',並且與mydata.Gest沒有關係。這就是我認爲它與'nls'有關的原因。 – MikeZ 2012-08-03 14:14:19