Skip to content

Commit

Permalink
BUG: brunnermunzel: fixed bugs pointed out in PR #7352
Browse files Browse the repository at this point in the history
  • Loading branch information
suzuki-shm authored and pv committed Apr 12, 2018
1 parent 82286a5 commit 2ad4b22
Showing 1 changed file with 10 additions and 9 deletions.
19 changes: 10 additions & 9 deletions scipy/stats/stats.py
Expand Up @@ -5267,27 +5267,28 @@ def brunnermunzel(x, y, alternative="two-sided", distribution="t", nan_policy='p
rankx_mean = np.mean(rankx)
ranky_mean = np.mean(ranky)

Sx = np.sum(np.power(rankcx - rankx - rankcx_mean + rankx_mean, 2))
Sx = np.sum(np.power(rankcx - rankx - rankcx_mean + rankx_mean, 2.0))
Sx /= nx - 1
Sy = np.sum(np.power(rankcy - ranky - rankcy_mean + ranky_mean, 2))
Sy = np.sum(np.power(rankcy - ranky - rankcy_mean + ranky_mean, 2.0))
Sy /= ny - 1

sigmax = Sx / np.power(nc - nx, 2)
sigmay = Sx / np.power(nc - ny, 2)
sigmac = nc * (sigmax / nx + sigmay / ny)
sigmax = Sx / np.power(nc - nx, 2.0)
sigmay = Sx / np.power(nc - ny, 2.0)

wbfn = nx * ny * (rankcy_mean - rankcx_mean)
wbfn /= (nx + ny) * np.sqrt(nx * Sx + ny * Sy)

if distribution == "t":
df = np.power(nx * Sx + ny * Sy,2)
df /= np.power(nx * Sx, 2) / (nx - 1) + np.power(ny * Sy,2) / (ny - 1)
df_numer = np.power(nx * Sx + ny * Sy, 2.0)
df_denom = np.power(nx * Sx, 2.0) / (nx - 1)
df_denom += np.power(ny * Sy, 2.0) / (ny - 1)
df = df_numer / df_denom
p = distributions.t.cdf(wbfn, df)
elif distribution == "norm":
elif distribution == "normal":
p = distributions.norm.cdf(wbfn)
else:
raise ValueError(
"distribution should be 't' or 'norm'")
"distribution should be 't' or 'normal'")

if alternative == "greater":
p = p
Expand Down

0 comments on commit 2ad4b22

Please sign in to comment.