Graviness Blog

算数・数学・科学・電脳・雑記・アホの順の密度で記事が構成されます.
<< April 2018 | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 >> ブログランキング・にほんブログ村へ
 
RECOMMEND
ビッグバン宇宙論 (上)
ビッグバン宇宙論 (上) (JUGEMレビュー »)
サイモン・シン, 青木 薫
RECENT COMMENT
  • 豊臣秀吉と曾呂利新左衛門から学ぶ数列の和
    優乃 (07/12)
  • 【誰か解いて】漸化式 a_(n+1) = f(n) * a_n ^ g(n) + h(n) の一般項
    優乃 (02/18)
  • 【誰か解いて】漸化式 a_(n+1) = f(n) * a_n ^ g(n) + h(n) の一般項
    S.S.+ (02/16)
  • 豊臣秀吉と曾呂利新左衛門から学ぶ数列の和
    坂井昭 (03/19)
  • d/dx(x↑↑n): 高さが定数のテトレーションの微分 - 数学的帰納法を用いる方法
    (09/30)
  • 全ての三角形は二等辺三角形
    優乃 (09/28)
  • 全ての三角形は二等辺三角形
    亀レス (09/28)
  • 全ての三角形は二等辺三角形
    優乃 (09/24)
  • 全ての三角形は二等辺三角形
    亀レス (09/23)
  • 【未解決】新しい演算子を創る
    $_ (09/10)
RECENT TRACKBACK
MOBILE
qrcode
PROFILE
無料ブログ作成サービス JUGEM
 
Brunner-Munzel検定のPython実装

Brunner-Munzel検定をPython(NumPy/SciPy使用)の実装です。

 

下記bmtest関数に2群のデータを入力すると、検定統計量W、p値(両側)、自由度、および帰無仮説で示される確率の推定値を出力します。

片側検定にする場合(どちらにずれているかが重要な場合)、出力されたp値を0.5倍して、Wの値が正負で判断します。

注意: 少しもかぶりもしない2群を入力すると、分散が0となって、検定統計量Wや自由度の計算中に Runtimewarning; divide by zero が出ます。

 

ref. Brunner-Munzel検定@奥村 晴彦さん

 

# coding: utf-8
#
# bmtest(x1, x2)
#   Calculate Brunner-Munzel-test scores.
#
#   Parameters:
#     x1, x2: array_like
#       Numeric data values from sample 1, 2.
#
#   Returns:
#     w:
#       Calculated test statistic.
#     p_value:
#       Two-tailed p-value of test.
#     dof:
#       Degree of freedom.
#     p:
#       "P(x1 < x2) + 0.5 P(x1 = x2)" estimates.
#
#   References:
#     * https://oku.edu.mie-u.ac.jp/~okumura/stat/brunner-munzel.html

import numpy
import scipy.stats

def bmtest(x1, x2):
    if isinstance(x1, numpy.ndarray): x1 = x1.tolist()
    if isinstance(x2, numpy.ndarray): x2 = x2.tolist()
    n1, n2 = len(x1), len(x2)
    R = scipy.stats.rankdata(x1 + x2, method = 'average')
    R1, R2 = R[:n1], R[n1:]
    r1_mean, r2_mean = numpy.mean(R1), numpy.mean(R2)
    p = (r2_mean - r1_mean) / (n1 + n2) + 0.5
    Ri1 = scipy.stats.rankdata(x1, method = 'average')
    Ri2 = scipy.stats.rankdata(x2, method = 'average')
    var1 = numpy.var([r - ri for r, ri in zip(R1, Ri1)], ddof = 1)
    var2 = numpy.var([r - ri for r, ri in zip(R2, Ri2)], ddof = 1)
    w = ((n1 * n2) * (r2_mean - r1_mean)) / ((n1 + n2) * numpy.sqrt(n1 * var1 + n2 * var2))
    dof = (n1 * var1 + n2 * var2) ** 2 / ((n1 * var1) ** 2 / (n1 - 1) + (n2 * var2) ** 2 / (n2 - 1))
    c = scipy.stats.t.cdf(abs(w), dof) if not numpy.isinf(w) else 0.0
    p_value = min(c, 1.0 - c) * 2.0
    return (w, p_value, dof, p)

 

JUGEMテーマ:コンピュータ

d/dx(x↑↑n): 高さが定数のテトレーションの微分

過去の記事「d/dx(x↑↑n): 高さが定数のテトレーションの微分 - 数学的帰納法を用いる方法」について、高さが定数のテトレーションの導関数を数学的帰納法を用いずに求めます。表記として、nx := x↑↑n を用います。

 

---

続きを読む >>
【誰か解いて】漸化式 a_(n+1) = f(n) * a_n ^ g(n) + h(n) の一般項

過去の二つの記事の漸化式を包含する a_(n+1) = f(n) * a_n ^ g(n) + h(n) の形式の一般項を求めたいです。

* 漸化式 a_(n+1) = f(n) * a_n + g(n) の一般項

* 漸化式 a_(n+1) = f(n) * a_n ^ g(n) の一般項

 

通勤時などで数ヶ月考えてますが、これがどうやっても解けない。。。恐らく初等関数の範囲では解けないのではないかと思います。

 

関連して、隣接三項間漸化式 a_(n+2) = f(n) * a_(n+1) + g(n) * a_n は、両辺を a_(n+1) で割ると、a_(n+2)/a_(n+1) = f(n) + g(n) * a_n/a_(n+1) となり、b_n = a_(n+1)/a_n とおけば、b_(n+1) = g(n) b_n ^ (-1) + f(n) となるため、本記事の漸化式が解ければ隣接三項間漸化式も解くことができます。

 

JUGEMテーマ:学問・学校

 

漸化式 a_(n+1) = f(n) * a_n ^ g(n) の一般項

あけましておめでとうございます。京都⇔宮崎の実家の移動中の数学の成果wを記事にします。

 

結論から記載すると、漸化式

a_(n+1) = f(n) * a_n ^ g(n)

の一般項が

a_n = {a_1^Π[i=1,n]g(i) * Π[j=1,n-1]{f(j)^Π[i=j+1,n]g(i)}}^{1/g(n)}

となることを示します。

 

---

続きを読む >>
階差数列から元の数列の一般項を求める問題では、n=1の場合分けが必要になることについて

高校数学で出題される次のような問題。

 

初項a1 = 1から始まる次の数列の一般項anを求めよ。
1, 2, 4, 7, 11, 16, 22, 29, 37

 

これを解く過程で場合分けが登場する。

 

各項の差をとると1, 2, 3, 4, 5, 6, 7, 8となり、階差数列の一般項はnと書ける。

 

ゆえに、n ≧ 2の場合

an = a1 + Σ[i=1,n-1]n

これを解いて

an = n(n-1)/2 + 1

ここで右辺にn = 1を代入すると

右辺 = 1 = a1 = 左辺

となり、n = 1のときも成り立つ。

 

最後のn = 1を確認するところで、成り立たない問題を見たことがない。だから、場合分けは必要なの?って話。以下で真面目に考えてみる。

続きを読む >>

(C) 2018 ブログ JUGEM Some Rights Reserved.