a_(n+2)=f(n) a_(n+1)+g(n) a_(n) の一般項
2019.08.15 Thursday
2019.04.30 Tuesday
確率密度関数f(x)=2xに従う確率変数の和と平均の確率密度関数の導出
興味のある分野等で統計学を活用してみると分かりますが、バシッと合うパターンが記載されていることは稀で、すぐに応用するための知識が必要になります。大学講義でも使用されるような書籍は濃過ぎだし、実用的か？というと非効率です。「プログラミングのための確率統計 」は応用し始めた人には丁度良いです。“プログラミングのための”となっていますが、プログラムは登場しません。“プログラミングのための”が同著の実用性の高さと知るべきに応える密度の濃さに一役買っています。また、感覚的にでも数式レベルで分からないと気持ちが悪いという人も必見です。本記事でも活用させて頂きました。 |
本記事では「中心極限定理の具体例 | 高校物理の備忘録」にある“具体例1”：の
が従う相対度数分布グラフの厳密な曲線、即ち
が従う確率密度関数を導出します。確率変数の変数変換、合成積（畳み込み）などを使用します。
JUGEMテーマ：学問・学校
2018.08.15 Wednesday
点の重なりを考慮した最適な点のサイズの導出
データを散布図等で表示するとき、データ数が多くて画面が点で塗り潰れてしまう経験はないでしょうか。
下図は10000個のデータを散布図で表示した例です。下図（左）は中央の帯の部分が点で塗り潰れています。描画領域を大きくするか、下図（右）のように点のサイズを小さく調節することで塗り潰れがなくなり、精確なデータ分析、そして美しく表示できます。
しかし、固定的なデータ数を扱う場合は問題となりませんが、データ数が変わるようなシチュエーションの場合はその都度調節する必要があり面倒です。
本記事では、データ数に応じて最適な点のサイズを自動計算する式とその導出を示します。結論を示すと、x個の点を描画するときの個々の点のサイズs(x)は次式で得られます。
s(x) = S1 * S∞ / (S1 * x + S∞)
S1: 点を1個描画するときの（x=1のときの）、1個の点による描画面積。
S∞: 総描画面積の上限。e.g. 画面の半分程度が埋まるのであればその面積。
例えば、S1=36px、S∞=3600pxの場合、s(x) = 3600 / (x + 100) となります。100点描画するとき個々の点のサイズはs(100)=18pxで描画され、同様にして1000点描画するときs(1000)=3.3px、10000点描画するときs(10000)=0.4pxとなります。
Pythonプログラムの例を示します。データ数len(data)に応じて点の描画サイズsを自動計算しています。
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
data = np.random.normal(size=1000) # normal dist.
S_1, S_m = 36, 3600
s = S_1 * S_m / (S_1 * len(data) + S_m) # see
plt.title(f"count={len(data)}")
plt.ylim((-4.0, 4.0))
plt.scatter(x=np.arange(len(data)), y=data, s=s)
plt.savefig(f"sample.{len(data)}.png")
plt.close()
以降で結論に至るまでの導出を示します。
JUGEMテーマ：コンピュータ
2018.05.05 Saturday
Smirnov-Grubbs検定を用いる外れ値除去のPython実装
スミルノフ−グラブス検定を用いる外れ値除去のPython（NumPy/SciPy使用）の実装です。
下記smirnov_grubbs関数にサンプルデータと有意水準α（両側100α％）を入力すると、外れ値を除去したサンプルデータと、抽出した外れ値データを出力します。
サンプルデータのうち平均値からの偏差が最も大きいデータを、有意水準（両側）でスミルノフ−グラブス検定し、帰無仮説を棄却した場合に当該データを外れ値とします。当該データを除いたサンプルデータに対し、外れ値がなくなるまでこの操作を繰り返します。
# coding: utf-8
#
# smirnov_grubbs(data, alpha)
# Generate outlier-removed data with Smirnov-Grubbs.
#
# Parameters:
# data: array_like
# Numeric data values.
# alpha: float
# Significance level. (two-sided)
#
# Returns:
# outlier-removed-data
# outlier-data
#
# e.g.
# data = [1, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 1000]
# alpha = 0.01
# smirnov_grubbs(data, alpha)
# ==> ([100, 101, 102, 103, 104, 105, 106, 107, 108, 109], [1000, 1])
import numpy as np
import scipy.stats as stats
def smirnov_grubbs(data, alpha):
x, o = list(data), []
while True:
n = len(x)
t = stats.t.isf(q=(alpha / n) / 2, df=n - 2)
tau = (n - 1) * t / np.sqrt(n * (n - 2) + n * t * t)
i_min, i_max = np.argmin(x), np.argmax(x)
myu, std = np.mean(x), np.std(x, ddof=1)
i_far = i_max if np.abs(x[i_max] - myu) > np.abs(x[i_min] - myu) else i_min
tau_far = np.abs((x[i_far] - myu) / std)
if tau_far < tau: break
o.append(x.pop(i_far))
return (np.array(x), np.array(o))
JUGEMテーマ：コンピュータ
2018.04.15 Sunday
Brunner-Munzel検定のPython実装
Brunner-Munzel検定のPython（NumPy/SciPy使用）の実装です。
下記bmtest関数に2群のデータを入力すると、検定統計量W、p値（両側）、自由度、および帰無仮説で示される確率の推定値を出力します。
片側検定にする場合（どちらにずれているかが重要な場合）、出力されたp値を0.5倍して、Wの値が正負で判断します。
注意: 少しもかぶりもしない2群を入力すると、分散が0となって、検定統計量Wや自由度の計算中に Runtimewarning; divide by zero が出ます。
# 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 as np
import scipy.stats as stats
def bmtest(x1, x2):
n1, n2 = len(x1), len(x2)
R = stats.rankdata(list(x1) + list(x2))
R1, R2 = R[:n1], R[n1:]
r1_mean, r2_mean = np.mean(R1), np.mean(R2)
Ri1, Ri2 = stats.rankdata(x1), stats.rankdata(x2)
var1 = np.var([r - ri for r, ri in zip(R1, Ri1)], ddof=1)
var2 = np.var([r - ri for r, ri in zip(R2, Ri2)], ddof=1)
w = ((n1 * n2) * (r2_mean - r1_mean)) / ((n1 + n2) * np.sqrt(n1 * var1 + n2 * var2))
dof = (n1 * var1 + n2 * var2) ** 2 / ((n1 * var1) ** 2 / (n1 - 1) + (n2 * var2) ** 2 / (n2 - 1))
c = stats.t.cdf(abs(w), dof) if not np.isinf(w) else 0.0
p_value = min(c, 1.0 - c) * 2.0
p = (r2_mean - r1_mean) / (n1 + n2) + 0.5
return (w, p_value, dof, p)
JUGEMテーマ：コンピュータ
2017.01.04 Wednesday
d/dx(x↑↑n): 高さが定数のテトレーションの微分
過去の記事「d/dx(x↑↑n): 高さが定数のテトレーションの微分 - 数学的帰納法を用いる方法」について、高さが定数のテトレーションの導関数を数学的帰納法を用いずに求めます。表記として、nx := x↑↑n を用います。
---
2017.01.04 Wednesday
【誰か解いて】漸化式 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テーマ：学問・学校
2016.11.20 Sunday
階差数列から元の数列の一般項を求める問題では、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) 2019 ブログ JUGEM Some Rights Reserved.
