2012-06-12 37 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 = "")