R中的嵌套for循环-最终结果的问题答案

【问题标题】：Nested for loops in R - Issue with final resultR中的嵌套for循环-最终结果的问题
【发布时间】：2012-06-14 23:14:27
【问题描述】：

当仅知道分布的矩时，我正在解决重建（或恢复）概率分布函数的问题。我已经用 R 编写了代码，虽然逻辑对我来说似乎是正确的，但我没有得到我想要的输出。

我试图用作近似（或重建或恢复）CDF 的方程式就是您在下图中看到的。我正在为等式的右侧编写代码，并将其等同于我在代码中称为 F 的向量。

可以在此处找到包含原始方程式的论文链接。

http://www.sciencedirect.com/science/article/pii/S0167715208000187

在论文中标记为等式（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 = "")

任何帮助都会在这里得到赞赏，因为使用我的代码获得的近似值和图形与原始曲线大不相同。

【问题讨论】：

仅供参考，我删除了块引用替代方案的内联代码降价。这可以通过将代码块缩进四个空格来完成，或者在突出显示要缩进的代码后点击{} 按钮。你能仔细检查一下我没有破坏代码的最初意图吗？
几乎不可能为您提供帮助，因为您的链接不可公开访问。因此，我既没有你的方程式，也没有你想要的结果。
当我运行您的代码时，我收到一个错误，即找不到 y。我猜 y 是一个向量，因为否则不是 k = 1 的长度？
看起来您正在尝试对阶乘使用自己的速记。我建议按照等式 2 的形式将它们写出来，直到你让代码工作为止。我发现这可以帮助我避免放错位置的括号。并不是说你有一个放错位置的括号。只是想提供一个有用的建议。另外，考虑在循环之外定义 mu，只是为了让代码更容易理解。我想我想弄清楚 9 是从哪里来的。

标签： r loops for-loop nested

【解决方案1】：

    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 = "")

【讨论】：

【解决方案2】：

您的问题似乎很明显。

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

你试图在 x 中得到不同的东西，是吗？那么，如果等式右边的 x 没有变化，而左边有一个使用非整数索引的赋值，你怎么能得到呢？

【讨论】：