2017-10-31 8 views
0

私はいくつかのデータのマンケンダルテストを実行しようとしています。私は次のリンク(https://github.com/mps9506/Mann-Kendall-Trend/blob/master/mk_test.py)のコードを使用して、結果が配列になるように少し変更し、p値(p)とτ値(z)のみを返します。Pythonで統計テストのために変な値を取得する

def mk_test(x, alpha=0.05): 
    """ 
    This function is derived from code originally posted by Sat Kumar Tomer 
    ([email protected]) 
    See also: http://vsp.pnnl.gov/help/Vsample/Design_Trend_Mann_Kendall.htm 
    The purpose of the Mann-Kendall (MK) test (Mann 1945, Kendall 1975, Gilbert 
    1987) is to statistically assess if there is a monotonic upward or downward 
    trend of the variable of interest over time. A monotonic upward (downward) 
    trend means that the variable consistently increases (decreases) through 
    time, but the trend may or may not be linear. The MK test can be used in 
    place of a parametric linear regression analysis, which can be used to test 
    if the slope of the estimated linear regression line is different from 
    zero. The regression analysis requires that the residuals from the fitted 
    regression line be normally distributed; an assumption not required by the 
    MK test, that is, the MK test is a non-parametric (distribution-free) test. 
    Hirsch, Slack and Smith (1982, page 107) indicate that the MK test is best 
    viewed as an exploratory analysis and is most appropriately used to 
    identify stations where changes are significant or of large magnitude and 
    to quantify these findings. 
    Input: 
     x: a vector of data 
     alpha: significance level (0.05 default) 
    Output: 
     trend: tells the trend (increasing, decreasing or no trend) 
     h: True (if trend is present) or False (if trend is absence) 
     p: p value of the significance test 
     z: normalized test statistics` 

    Examples 
    -------- 
     >>> x = np.random.rand(100) 
     >>> trend,h,p,z = mk_test(x,0.05) 
    """ 
    n = len(x) 

    # calculate S 
    s = 0 
    for k in range(n-1): 
     for j in range(k+1, n): 
      s += np.sign(x[j] - x[k]) 

    # calculate the unique data 
    unique_x = np.unique(x) 
    g = len(unique_x) 

    # calculate the var(s) 
    if n == g: # there is no tie 
     var_s = (n*(n-1)*(2*n+5))/18 
    else: # there are some ties in data 
     tp = np.zeros(unique_x.shape) 
     for i in range(len(unique_x)): 
      tp[i] = sum(x == unique_x[i]) 
     var_s = (n*(n-1)*(2*n+5) - np.sum(tp*(tp-1)*(2*tp+5)))/18 

    if s > 0: 
     z = (s - 1)/np.sqrt(var_s) 
     #result = (s - 1)/np.sqrt(var_s) 
    elif s == 0: 
     z = 0 
     #result = 0 
    elif s < 0: 
     z = (s + 1)/np.sqrt(var_s) 
     #result = (s + 1)/np.sqrt(var_s) 

    # calculate the p_value 
    p = 2*(1-norm.cdf(abs(z))) # two tail test 
    result= np.append(p,z) 
    h = abs(z) > norm.ppf(1-alpha/2) 

    return np.array(result) 

次に、次のコードを使用してテストを行っています。

out = np.empty((0)) 
for i in range(145): 
    for j in range(192): 
     out1 = mk_test(yrmax[:,i,j], alpha=0.05) 
     out = np.append(out, out1, axis=0) 

は、私は-1と1の間のzの値を取得することを期待するので、テストの実行中に何かがどこか間違っていると思いますが、私は1より大きいいくつかの値を取り戻すてることは起こって符号化です間違っているか、私はzが何であるか誤解しています。それは実際にタウではなく、なぜ私が期待していない価値を得ているのですか?

+1

この質問はあなたの統計のcalcuation /実装についてであるならば、あなたはhttps://stats.stackexchange.com/ –

+0

でより良い運を持っているかもしれない、それは間違っていると何か起こっていたかどうかはわかりませんでしたコーディングかどうか。私もそこに投稿します、ありがとう! –

答えて

関連する問題