2017-02-28 40 views
1

我想計算安德森達令測試發現here。我遵循維基百科上的步驟,並確保在計算我正在測試的數據的平均值和標準偏差時,使用MATLAB表示爲X。另外,我使用了一個叫phi的函數來計算標準正常CDF,我也測試了這個函數以確保它是正確的。現在,當我實際計算A平方(在維基百科中表示,我在C++中將其表示爲A)時,似乎出現了問題。安德森達林測試C++

這裏是我的功能,我對安德森 - 達林檢驗進行:

void Anderson_Darling(int n, double X[]){ 
    sort(X,X + n); 
    // Find the mean of X 
    double X_avg = 0.0; 
    double sum = 0.0; 
    for(int i = 0; i < n; i++){ 
     sum += X[i]; 
    } 
    X_avg = ((double)sum)/n; 


    // Find the variance of X 
    double X_sig = 0.0; 
    for(int i = 0; i < n; i++){ 
     X_sig += (X[i] - X_avg)*(X[i] - X_avg); 
    } 
    X_sig /= n; 


    // The values X_i are standardized to create new values Y_i 
    double Y[n]; 
    for(int i = 0; i < n; i++){ 
     Y[i] = (X[i] - X_avg)/(sqrt(X_sig)); 
     //cout << Y[i] << endl; 
    } 

    // With a standard normal CDF, we calculate the Anderson_Darling Statistic 
    double A = 0.0; 
    for(int i = 0; i < n; i++){ 
     A += -n - 1/n *(2*(i) - 1)*(log(phi(Y[i])) + log(1 - phi(Y[n+1 - i]))); 
    } 
    cout << A << endl; 
} 

注意,我知道,安德森 - 達林式(A方)與i = 1開始i = n,雖然當我改變了索引使它在C++中工作,我仍然得到相同的結果而不改變索引。

我在C++中獲得的價值是:

-4e+006 

我應該得到的值,在MATLAB收到的是:

0.2330 

任何建議都非常讚賞。

這裏是我的全部代碼:

#include <iostream> 
#include <math.h> 
#include <cmath> 
#include <random> 
#include <algorithm> 
#include <chrono> 

using namespace std; 

double *Box_Muller(int n, double u[]); 
double *Beasley_Springer_Moro(int n, double u[]); 
void Anderson_Darling(int n, double X[]); 
double phi(double x); 

int main(){ 
    int n = 2000; 
    double Mersenne[n]; 
    random_device rd; 
    mt19937 e2(1); 
    uniform_real_distribution<double> dist(0, 1); 
    for(int i = 0; i < n; i++){ 
     Mersenne[i] = dist(e2); 
    } 

    // Print Anderson Statistic for Mersenne 6a 
    double *result = new double[n]; 
    result = Box_Muller(n,Mersenne); 
    Anderson_Darling(n,result); 




    return 0; 
} 

double *Box_Muller(int n, double u[]){ 
    double *X = new double[n]; 
    double Y[n]; 
    double R_2[n]; 
    double theta[n]; 
    for(int i = 0; i < n; i++){ 
     R_2[i] = -2.0*log(u[i]); 
     theta[i] = 2.0*M_PI*u[i+1]; 
    } 
    for(int i = 0; i < n; i++){ 
     X[i] = sqrt(-2.0*log(u[i]))*cos(2.0*M_PI*u[i+1]); 
     Y[i] = sqrt(-2.0*log(u[i]))*sin(2.0*M_PI*u[i+1]); 
    } 
    return X; 
} 

double *Beasley_Springer_Moro(int n, double u[]){ 
    double y[n]; 
    double r[n+1]; 
    double *x = new double(n); 
    // Constants needed for algo 
    double a_0 = 2.50662823884;  double b_0 = -8.47351093090; 
    double a_1 = -18.61500062529; double b_1 = 23.08336743743; 
    double a_2 = 41.39119773534; double b_2 = -21.06224101826; 
    double a_3 = -25.44106049637; double b_3 = 3.13082909833; 

    double c_0 = 0.3374754822726147; double c_5 = 0.0003951896511919; 
    double c_1 = 0.9761690190917186; double c_6 = 0.0000321767881768; 
    double c_2 = 0.1607979714918209; double c_7 = 0.0000002888167364; 
    double c_3 = 0.0276438810333863; double c_8 = 0.0000003960315187; 
    double c_4 = 0.0038405729373609; 

    // Set r and x to empty for now 
    for(int i = 0; i <= n; i++){ 
     r[i] = 0.0; 
     x[i] = 0.0; 
    } 
    for(int i = 1; i <= n; i++){ 
     y[i] = u[i] - 0.5; 
     if(fabs(y[i]) < 0.42){ 
      r[i] = pow(y[i],2.0); 
      x[i] = y[i]*(((a_3*r[i] + a_2)*r[i] + a_1)*r[i] + a_0)/((((b_3*r[i] + b_2)*r[i] + b_1)*r[i] + b_0)*r[i] + 1); 
     }else{ 
      r[i] = u[i]; 
      if(y[i] > 0.0){ 
       r[i] = 1.0 - u[i]; 
       r[i] = log(-log(r[i])); 
       x[i] = c_0 + r[i]*(c_1 + r[i]*(c_2 + r[i]*(c_3 + r[i]*(c_4 + r[i]*(c_5 + r[i]*(c_6 + r[i]*(c_7 + r[i]*c_8))))))); 
      } 
      if(y[i] < 0){ 
       x[i] = -x[i]; 
      } 
     } 
    } 
    return x; 
} 

    double phi(double x){ 
    return 0.5 * erfc(-x * M_SQRT1_2); 
} 


void Anderson_Darling(int n, double X[]){ 
    sort(X,X + n); 
    // Find the mean of X 
    double X_avg = 0.0; 
    double sum = 0.0; 
    for(int i = 0; i < n; i++){ 
     sum += X[i]; 
    } 
    X_avg = ((double)sum)/n; 


    // Find the variance of X 
    double X_sig = 0.0; 
    for(int i = 0; i < n; i++){ 
     X_sig += (X[i] - X_avg)*(X[i] - X_avg); 
    } 
    X_sig /= (n-1); 


    // The values X_i are standardized to create new values Y_i 
    double Y[n]; 
    for(int i = 0; i < n; i++){ 
     Y[i] = (X[i] - X_avg)/(sqrt(X_sig)); 
     //cout << Y[i] << endl; 
    } 

    // With a standard normal CDF, we calculate the Anderson_Darling Statistic 
    double A = -n; 
    for(int i = 0; i < n; i++){ 
     A += -1.0/(double)n *(2*(i+1) - 1)*(log(phi(Y[i])) + log(1 - phi(Y[n - i]))); 
    } 
    cout << A << endl; 
} 
+2

在調試器中瀏覽代碼並找出發生了什麼。 –

+0

@AnonMail好的我正在使用Eclipse,通過確保我具有用於計算X的平均值和標準偏差的相同值進行調試。您是否建議使用Eclipse調試器? – Scooby

+0

在附註中,你爲什麼要重寫自己的'Phi'?我在過去做過,但後來轉而使用'std :: erfc'和'std :: erf',這非常準確和快速。 –

回答

3

讓我猜,你的n是2000.對嗎? 這裏的主要問題是當你在最後一個表達式中做1/n。 1是一個int並且ao是n。當你用n除1時,它執行整數除法。現在1除以任意數字> 1在整數除法下爲0(認爲如果它只保留商的整數部分,那麼你需要做的是通過寫入1 /(double)n來將n寫成雙精度數n

其餘的就應該可以正常

摘要從討論 -

  1. 索引到Y []應該是我和n-1-I分別爲
  2. n應不會在循環,但僅添加一次。 。

小修正像計算方差時將除數改爲n而不是n-1。

+0

我看,我試過這樣做,但我得到這個值-3.998e + 006和Matlab得到0.2330,所以我不知道什麼是錯的 – Scooby

+0

得到了你的錯誤。這不是一個編碼錯誤,但是您在翻譯公式時犯了一個錯誤。 A是 - n - S。 S是總和。你每次都添加了-n。結果你獲得了巨大的負面價值。初始化A = -n並將其從循環中移除。 –

+0

對不起,你可以更具體地說明我需要做什麼,我完全不理解你的評論 – Scooby

2

你有整數除法這裏:

A += -n - 1/n *(2*(i) - 1)*(log(phi(Y[i])) + log(1 - phi(Y[n+1 - i]))); 
      ^^^ 

1/n爲零時n > 1 - 你需要這種改變,如:1.0/n

A += -n - 1.0/n *(2*(i) - 1)*(log(phi(Y[i])) + log(1 - phi(Y[n+1 - i]))); 
      ^^^^^ 
+0

謝謝,我這樣做了,但是我得到了-3.998e + 006而不是正確的值。由於A平方的公式是從i = 1開始的,我從i = 0開始,因爲這就是數組如何在C++中編入索引,是否需要更改A的公式?也許讓我所有的i + 1? – Scooby

+0

我對你使用的公式不熟悉,但是從你所說的聽起來你需要修復索引,即使用'(i + 1)'而不是'i'。 –

+0

這很奇怪,我仍然得到相同的值-3.998e + 006,當它應該是.2330。 – Scooby