2

確率分布関数の再構築 (または回復) の問題を解決している最中ですが、分布の瞬間のみがわかっています。私は R でコードを書きましたが、ロジックは正しいように思えますが、必要な出力が得られません。

近似 (または再構築または回復) 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 = "")

私のコードを使用して得られた近似とグラフは元の曲線とは大きく異なるため、ここで助けていただければ幸いです。

4

2 に答える 2

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 の変化がなく、左辺に整数以外のインデックスを使用した割り当てがある場合、どうすればそれを得ることができるでしょうか?

于 2012-06-12T08:32:59.390 に答える
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 = "")
于 2012-06-12T11:31:50.787 に答える