pythonで完全数と友愛数と婚約数を求める

以前の記事の続きとして、せっかくなので友愛数・完全数・婚約数を求めてみたい。

約数関数（再掲）

divisor(n)として約数関数を使うので、再掲しておく。
sympyを使う。

import sympy
def divisor(n):
    factors = sympy.factorint(n)
    rst = 1
    for i,j in factors.iteritems():
        rst *= (pow(i,j+1) - 1)/(i-1)
    return rst

完全数

完全数（かんぜんすう，英: perfect number）とは、自分自身を除く正の約数の和に等しくなる自然数のことである。完全数の最初の3個は 6 (= 1 + 2 + 3)、28 (= 1 + 2 + 4 + 7 + 14)、496 (= 1 + 2 + 4 + 8 + 16 + 31 + 62 + 124 + 248) である。「完全数」は「万物は数なり」と考えたピタゴラスが名付けた数の一つであることに由来する[1]が、彼がなぜ「完全」と考えたのかについては何も書き残されていないようである
wikipedia:完全数

max = 100000
rst = []
for n in range(1,max+1):
    if divisor(n) == 2*n:
        rst.append(n)

• 結果

rst = 
[6, 28, 496, 8128]

友愛数

友愛数（ゆうあいすう、英: amicable numbers）とは、異なる 2 つの自然数の組で、自分自身を除いた約数の和が、互いに他方と等しくなるような数をいう。親和数（しんわすう）とも呼ばれる。
最小の友愛数の組は (220, 284) である。
220 の自分自身を除いた約数は、1, 2, 4, 5, 10, 11, 20, 22, 44, 55, 110 で、和は 284 となる。一方、284 の自分自身を除いた約数は、1, 2, 4, 71, 142 で、和は 220 である。
wikipedia:友愛数

max = 100000
rst = []
for n in range(1,max+1):
    pair = divisor(n) - n
    if n < pair and pair <= max:
        if divisor(pair) - pair == n:
            rst.append([n,pair])

• 結果

rst = 
[[220, 284],
 [1184, 1210],
 [2620, 2924],
 [5020, 5564],
 [6232, 6368],
 [10744, 10856],
 [12285, 14595],
 [17296, 18416],
 [63020, 76084],
 [66928, 66992],
 [67095, 71145],
 [69615, 87633],
 [79750, 88730]]

婚約数

婚約数（こんやくすう、英: betrothed numbers）とは、異なる 2 つの自然数の組で、1 と自分自身を除いた約数の和が、互いに他方と等しくなるような数をいう。準友愛数（じゅんゆうあいすう、quasi-amicable numbers）とも呼ばれる。
最小の婚約数の組は (48, 75) である。
48 の 1 と自分自身を除いた約数は、2, 3, 4, 6, 8, 12, 16, 24 で、和は 75 となる。一方、75 の 1 と自分自身を除いた約数は、3, 5, 15, 25 で、和は 48 である。
wikipedia:婚約数

max = 100000
rst = []
for n in range(1,max+1):
    pair = divisor(n) - n - 1
    if n < pair and pair <= max:
        if divisor(pair) - pair - 1 == n:
            rst.append([n,pair])

• 結果

rst = 
[[48, 75],
 [140, 195],
 [1050, 1925],
 [1575, 1648],
 [2024, 2295],
 [5775, 6128],
 [8892, 16587],
 [9504, 20735],
 [62744, 75495]]

sympyで素因数分解ができるからオイラー関数と約数関数を書いてみた

pythonのライブラリの一つsympyを使えば、簡単に素因数分解ができるということを知った。

import sympy
sympy.factorint(1000)
#{2: 3, 5: 3}

ちなみに、因数分解も簡単にできる!!
凄い。

import sympy
x = sympy.Symbol('x')

eqn = x**2 - 3*x + 2
print(sympy.factor(eqn))
# (x - 2)*(x - 1)

これを使って、オイラー関数と約数関数を書いてみた。

オイラー関数

自然数に対して、未満でと互いに素な自然数の個数をで表し、オイラー関数と呼ぶ。
例えば、未満でと互いに素な自然数はなので

def euler(n):
    factors = sympy.factorint(n)
    rst = 1
    for i,j in factors.iteritems():
        rst *= (pow(i,j) - pow(i,j-1))
    return rst

約数関数

自然数の約数の和をと表し、約数関数と呼ぶ。
例えば、の約数はなので

def divisor(n):
    factors = sympy.factorint(n)
    rst = 1
    for i,j in factors.iteritems():
        rst *= (pow(i,j+1) - 1)/(i-1)
    return rst

コードの説明

コードだけ見てもちんぷんかんぷんだろうけれど、それぞれ関数のこれらの性質を利用している。

の素因数分解が次のように与えられているとき

オイラー関数、約数関数はそれぞれ次のように求まる

数え上げ

n	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19
	1	1	2	2	4	2	6	4	6	4	10	4	12	6	8	8	16	6	18
	1	3	4	7	6	12	8	15	13	18	12	28	14	24	24	31	18	39	20

参考：
wikipedia:オイラー関数
wikipedia:約数関数
オイラー関数に関して：
RSA暗号を実装してみる（理論編） - 技術メモ

関連記事：

pythonで一元配置分散分析（one way ANOVA)

一元配置分散分析とは

「3つ以上の群があった時に、果たしてそれらの群の平均は等しいと言えるかどうか。」
という検定。
集団の分布が正規性を持つことが前提となっている。*1
※すべての組み合わせペアでt検定を適用するのは間違いなので注意*2

基本的な発想は、
「集団間の分散と集団内のばらつきを比べて、集団間のばらつきの方が大きいと見なせれば集団間に違いがある。」
と考えるという仕組みになっている。

方法

帰無仮説と対立仮設は以下のようになる。

• 帰無仮説：集団間の平均に違いはない
• 対立仮設：少なくとも1つの集団の平均が他の集団と異なる

上で説明したような
「集団内の分散から見た集団間のばらつきの大きさ」
をF値と呼び、これがF分布に従うことを利用して検定を行う。

F値が有意に大きい見なされると帰無仮説が棄却され、「集団の平均がすべて同じとは言えない」と言える。

Pythonのコード

例によってscipyが使える。

scipy.stats.f_oneway(*args)

argsには集団のデータが入る。
p値はpvalueで取り出す。

テストコード

平均が同じ集団

まず、帰無仮説が採択される場合を見てみる。
例えば、下の3つを比較してみよう。

1. 平均10　分散5　標本数90
2. 平均10　分散10　標本数120
3. 平均10　分散20　標本数100

import matplotlib.pyplot as plt
import numpy as np

a1 = np.random.normal(10,5,size=90)
a2 = np.random.normal(10,10,size=120)
a3 = np.random.normal(10,20,size=100)
bins = np.arange(-30,50,5)
plt.hist(a1, bins=bins, alpha=0.3)
plt.hist(a2, bins=bins, alpha=0.3)
plt.hist(a3, bins=bins, alpha=0.3)
result = scipy.stats.f_oneway(a1,a2,a3)
print(result.pvalue)
plt.show()

0.352570030224

この場合、
有意水準5％とすると
0.352>0.05
なので、集団間の平均値に差異はないことが分かった。

平均が異なる集団

次に、帰無仮説が棄却される場合を見てみる。
上のデータを少し変えて

1. 平均10　分散5　標本数90
2. 平均15　分散10　標本数120
3. 平均10　分散20　標本数100

としてみる

import matplotlib.pyplot as plt
import numpy as np

a1 = np.random.normal(10,5,size=90)
a2 = np.random.normal(15,10,size=120)
a3 = np.random.normal(10,20,size=100)
bins = np.arange(-30,50,5)
plt.hist(a1, bins=bins, alpha=0.3)
plt.hist(a2, bins=bins, alpha=0.3)
plt.hist(a3, bins=bins, alpha=0.3)
result = scipy.stats.f_oneway(a1,a2,a3)
print(result.pvalue)
plt.show()

0.00362625959975

この場合は、
有意水準5％とすると
0.004<0.05より
帰無仮説は棄却される。
「少なくとも1つの集団の平均値が異なる」ことがわかる。

参考サイト
立教大学の授業スライド：わかりやすくまとまってある。リンク切れ
統計検定を理解せず使っている人のためにⅢ
scipy.stats.f_oneway — SciPy v0.19.0 Reference Guide

関連記事：

*1:ノンパラメトリックの場合はKruskal-Wallis検定を使うらしい。勉強中。scipy.stats.mstats.kruskalwallis — SciPy v0.14.0 Reference Guide

*2:https://www.numberanalytics.com/