Rのコードをpythonに書き換えたい

Question

Rのコードをpythonに書き換えたいと思っています。
実装したい内容は、以下で、ziは標準正規分布の累積分布関数をあらわします。
アンダーソンダーリング検定を1から実行して、最終的にA*>0.752 がシミュレーション回数のうち何個あるかで、検出力を出力したいと考えています。

Rのコードは以下です。

n = 50
sim<-100
count<-0

MAD<-function(){
    for(i in 1:sim){
    #rdata<-rlogis(n,0,1)
    rdata<-rt(n,1)#累積分布関数に使用する分布から乱数をとる。
    data<-sort(rdata)#順序統計量(小さい順に並べかえる)
    cdf<-pnorm(data, mean(data), sd(data))(b11f778f7ed7ac3a11db5dfa9fb94c49.png)
    S <- 0
    for(i in 1:n){
        S <- S + (2*i-1) * (log(cdf[i])+log(1 - cdf[n + 1 - i]))
    }
    AD <- ( - n - S/n)*(1 + 0.75/n + 2.25/n^2)
    #print(b1)# 統計量
    if(AD >= 0.752){
        count <- count + 1
    }
}
print(count/sim)#検出力
}
MAD()

自力で書き換えたpythonは以下です。しかし、print(MAD)の部分が1より大きい値になってしまい、どこが誤っているのかわからない状態です。

def MAD():
    count = 0
    sim = 100
    for s in range(1,sim):
        rdata = np.random.standard_t(1, size=50)
        data = sorted(rdata)
        n = len(data)
        mu = np.mean(data)
        seg = np.std(data)
        # cdfで累積分布関数を生成
        # F0　標準正規分布の累積分布関数
        S = 0
        for i in range(1,n):
            X = data[i]
            F0 = stats.norm.cdf(x=X, loc=mu, scale=seg)
            XX = data[n- i]
            FF = stats.norm.cdf(x=XX, loc=mu, scale=seg)
            S += (2 * i - 1)*(math.exp(F0) + math.exp(1-FF))
         MAD = (- n - S/n) *(1 + 0.75 / n + 2.25 / n**2)
         print(MAD)#統計量
         if MAD >= 0.752:
            count += 1
         else:
             pass
    power = count / sim
    print(count)
    print(power)

MAD()

Rを書き換えたpythonのコードを教えて頂きたいです。
よろしくお願い致します。

Comments

[KojiDoi]

Sの値は間違いなく算出出来てるのでしょうか。python版とR版で微妙に式の立て方が違うようですが。

[emiime]

R版が正しく、自分で考えたのは、python版なのですが、自分では同じだと思って、pythonを組んでみたのですが、どこが違っているのかがわからなく詰まってしまいました。
ただ、python版のXX = data[n - i] の部分は エラーが出てしまったため、一旦XX = data[n - i]としましたが、XX = data[n + 1 - i] に直そうと思っています。

Rと比べた時にpythonはどこが異なっているのかでもわかれば教えて頂きたいです。

[KojiDoi]

>どこが異なっているのか
煩雑なのでコードはきちんと追ってませんけど、どうみてもSの定義の式は違っていますよね。

なぜ変えているのですか? まずは、極力R版に寄せて組んでみるべきです。変数名もできるだけ同じにしましょう。そして、Sの値を出力して、両版で違いが出ないか見るのです。Sが問題なさそうなら、その前のcdfとやらを比較して違い長いかチェックする。そういう基本的なチェックがいい加減なまま悩み続けるのは時間の無駄です。

[bsdfan]

まずは乱数を使わずに、同じ配列を入れた時に同じ値になるかを、それぞれの行ごとにチェックしながらやるといいと思います。

[emiime]

乱数を使わずにやってみます！ありがとうございます！

Answer 1

ちゃんと中まで確認していませんが、Rの式では、2つめのループの外にADの計算がありますが、pythonのものは2つめのループの中に入っています。

2つめのループの中で(M)ADを計算して条件に合う数を数えているのですから、simより大きくなってしまうのでしょう。

Comments

[emiime]

for文の外に、MAD =の式を出してみたのですが、(上記編集の通り)そうすると、1つ目のfor文でSが変化するのに対して、MADの式に当てはまるSが変化しないことになってしまうと考え、
MADをループから出す方法？の書き方がよくわからない、と詰まってしまいました。
もし、追加のアドバイス頂けると嬉しいです。

[TakaiY]

僕も詳しいわけではありませんが、
ループの中にある、この式でSは積算されていますよね。

S <- S + (2*i-1) * (log(cdf[i])+log(1 - cdf[n + 1 - i]))

その積算結果を使っていますから、それでいいのではないかと思いますよ。

pythonの処理がRの処理と同等なのかどうかはわかりません。
ですが、pythonの式は引いているのがちょっと気にはなります。
S -= (2 * i - 1)*(math.exp(F0) + math.exp(1-FF))