2012-06-12 40 views
2

我正在解決重建(或恢復)概率分佈函數時只知道分佈的矩的問題。我已經爲它編寫了代碼,雖然邏輯看起來是正確的,但我沒有得到我想要的輸出。在R中嵌套for循環 - 問題與最終結果

我試圖用作近似(或重建或恢復)CDF的公式是您在下圖中看到的。我正在編寫等式右邊的代碼,並將其與我在代碼中稱爲F的向量等同起來。

包含原始公式的紙張鏈接可以在此處找到。

它被標記爲(2)在紙方程。

這裏是我寫的:

#R Codes: 
alpha <- 50 
T <- 1 
x <- seq(0, T, by = 0.1) 

# Original CDF equation 
Ft <- (1-log(x^3))*(x^3) 
plot(x, Ft, type = "l", ylab = "", xlab = "") 

# Approximated CDF equation using Moment type reconstruction 
k<- floor(alpha*y/T) 
for(i in 1:length(k)) 
{ 
for(j in k[i]:alpha) 
{ 
F[x+1] <- (factorial(alpha)/(factorial(alpha-j)*factorial(j-k)*factorial(k)))*(((-1)^(j-k))/(T^j))*((9/(j+3))^2) 
} 
} 
plot(x[1:7], F, type = "l", ylab = "", xlab = "") 

任何幫助將這裏讚賞,因爲使用我的代碼獲得的逼近和圖形是從原來的曲線非常不同的代碼。

+1

僅供參考,我刪除了您的內聯代碼減價塊的報價替代。這可以通過用四個空格縮進代碼塊來實現,或者在突出顯示要縮進的代碼後觸擊'{}'按鈕。你可以仔細檢查一下,我沒有破壞代碼的初始意圖嗎? – Chase

+2

這幾乎不可能幫助你,因爲你的鏈接不公開。因此,我既沒有你的等式,也沒有你想要的結果。 – Andrie

+0

當我運行你的代碼時,我得到一個錯誤,找不到y。我猜y是一個向量,因爲否則不是k = 1的長度? –

回答

1

看來很清楚你的問題在這裏。

F[x+1] <- (factorial(alpha)/(factorial(alpha-j)*factorial(j-k)*factorial(k)))*(((-1)^(j-k))/(T^j))*((9/(j+3))^2) 

你試圖讓東西變化在x,是嗎?所以,如果這個方程的右邊在x中沒有任何變化,而左邊有一個使用非整數索引的賦值,那麼怎麼能得到呢?

0
alpha <- 30 #In the exemple you try to reproduce, they use an alpha of 30 if i understood correctly (i'm a paleontologist not a mathematician so this paper's way beyond my area of expertise :)) 

    tau <- 1 #tau is your T (i changed it to avoid confusion with TRUE) 
    x <- seq(0, tau, by = 0.001) 
    f<-rep(0,length(x)) #This is your F (same reason as above for the change). 
    #It has to be created as a vector of 0 before your loop since the whole idea of the loop is that you want to proceed by incrementation. 

    #You want a value of f for each of your element of x so here is your first loop: 
    for(i in 1:length(x)){ 

    #Then you want the sum for all k going from 1 to alpha*x[i]/tau: 
     for(k in 1:floor(alpha*x[i]/tau)){ 

    #And inside that sum, the sum for all j going from k to alpha: 
      for(j in k:alpha){ 

    #This sum needs to be incremented (hence f[i] on both side) 
       f[i]<-f[i]+(factorial(alpha)/(factorial(alpha-j)*factorial(j-k)*factorial(k)))*(((-1)^(j-k))/(tau^j))*(9/(j+3)^2) 
       } 
      } 
     } 

    plot(x, f, type = "l", ylab = "", xlab = "")