Go to Qiita Advent Calendar 2024 Top
LoginSignup
5
4

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

数値計算Advent Calendar 2020

Day 11
@Snakeseed

統計で使える身近な関数をPythonで作ってみた

Last updated at Posted at 2020-12-18

##目次

  • この記事を書こうと思ったきっかけ
  • 実行環境
  • 平均
  • レンジ
  • 中央値1
  • 中央値2(昇順、降順の時のみ)
  • 標準偏差
  • 分散
  • 累積相対度数
  • 変動係数
  • 2変量解析
    • 共分散
    • 相関係数
  • おまけ編

##この記事を書こうと思ったきっかけ
図書館で統計学の本を読んでいたところ、「Pythonを使って(標準偏差などを)求められるのでは」と思ったこと。また、Pythonの実力をあげたかったので、気分転換も兼ねてこの記事を書きました。

##実行環境
macOS Big Sur
11.1
paiza.IO(Python3)で実行しています。

##平均
####プログラム

def mean(l):
    return sum(l)/len(l)

####使用例

print(mean([3,4,5,7,6]))
#5.0
print(mean([3,4,11,7,6]))
#6.2

####注意点

  • 戻り値はfloat型
  • 昇順でないリストにも使える点
  • 数がとてつもなく大きくなると下のように「e」で省略されてしまう点
print(mean([3,4,1111111111111111111111111111111,7,6]))
#2.222222222222222e+29

##レンジ
####プログラム

def Range(l):
    return max(l)-min(l)

昇順でなくても良い。

####注意点

  • 組込み変数のrangeと重ならないように関数名の1文字目を大文字にしています。

##中央値1
####プログラム

def chuuouti(l):
    l2 = sorted(l)
    if len(l2)%2:
        return l2[int(len(l)/2)]
    else:
        return (l2[int(len(l)/2)-1]+l2[int(len(l)/2)])/2

####使用例

print(chuuouti([1,4,2,3]))

####注意点

  • 要素数が偶数の時は真ん中の2つの要素の平均をとる。
  • 昇順でなくても良い(しかし、sortedの処理に少し時間がかかるので、昇順の時は「中央値2」を使うと良い。

####例題
第四回 アルゴリズム実技検定　A - 中央値

##中央値2(昇順、降順の時のみ)
####プログラム

def chuuouti(l):
    if len(l)%2:
        return l[int(len(l)/2)]
    else:
        return (l[int(len(l)/2)-1]+l[int(len(l)/2)])/2

####注意点

  • 昇順、降順の時以外は使えない
  • 長さが偶数である時は戻り値がfloat型だが、長さが奇数の時は戻り値の型が中央値と一致する点。(下)
#長さが偶数の時
print(chuuouti([3,2,2,1]))
#2.0
print(chuuouti([4,3,2,1]))
#2.5

#長さが奇数の時
print(chuuouti([1,3,7]))
#3
print(chuuouti([1,3.2,7]))
#3.2

##標準偏差
####プログラム

from math import sqrt
def deviation_value(l):
    mean = sum(l)/len(l)#平均
    deviation_square = [(l[i]-mean)**2 for i in range(len(l))]#偏差平方
    deviation_square_sum = sum(deviation_square)#偏差平方和
    return sqrt(deviation_square_sum/len(l))#偏差値

####使用例

print(deviation_value([10,15,5,-22,-18]))
#15.086417732516887

上の使用例では平均が{10+15+5+(-22)+(-18)}/5より-2.
よって、偏差平方は144,289,49,400,256.この和は1138.
これを要素数(5個)で割り、平方根をとれば、標準偏差になる。

##累積相対度数
####プログラム

def cumulative_relative_frequency(x):
    wa = sum(x)#全体の和
    l = [0]
    y = [x[i]/wa for i in range(len(x))]
    #ここから累積和
    for i in range(len(y)):
        l.append(l[-1]+y[i])
    del l[0]
    return l

####使用例

print(cumulative_relative_frequency([1,3,7,6,9]))
#[0.038461538461538464, 0.15384615384615385, 0.4230769230769231, 0.6538461538461539, 1.0]

##変動係数
####プログラム

from math import sqrt
def coefficient_of_variation(l):
    mean = sum(l)/len(l)#平均
    deviation_square = [(l[i]-mean)**2 for i in range(len(l))]#偏差平方
    deviation_square_sum = sum(deviation_square)#偏差平方和
    return sqrt(deviation_square_sum/len(l))/mean

####使用例

print(coefficient_of_variation([4,5,3,5]))
#https://bellcurve.jp/statistics/course/19515.htmlの例。

本当は0.2であるが、0.195と出力しているため、少し誤差があるとみられる。この誤差は、平方根(sqrt)の時に発生したものとみられる。

#2変量解析
##共分散
####プログラム

def bivariate_analysis(x,y):
    mean_x = sum(x)/len(x)
    mean_y = sum(y)/len(y)
    deviation_x = [x[i]-mean_x for i in range(len(x))]
    deviation_y = [y[i]-mean_y for i in range(len(y))]
    product_of_deviations = [deviation_x[i]*deviation_y[i] for i in range(len(x))]
    return sum(product_of_deviations)/len(x)

####使用例

print(bivariate_analysis([4,7],[3,6]))

この関数を使うことで、2つのリストの共分散が求められる。

####メリット

  • 正確な値を求められる。

##相関係数
####プログラム

def correlation_coefficient(x,y):
    return bivariate_analysis(x,y)/deviation_value(x)/deviation_value(y)

####標準偏差と共分散を入れたプログラム

from math import sqrt

def deviation_value(l):#標準偏差
    mean = sum(l)/len(l)#平均
    deviation_square = [(l[i]-mean)**2 for i in range(len(l))]#偏差平方
    deviation_square_sum = sum(deviation_square)#偏差平方和
    return sqrt(deviation_square_sum/len(l))#標準偏差
    
    
def bivariate_analysis(x,y):#共分散
    mean_x = sum(x)/len(x)
    mean_y = sum(y)/len(y)
    deviation_x = [x[i]-mean_x for i in range(len(x))]
    deviation_y = [y[i]-mean_y for i in range(len(y))]
    product_of_deviations = [deviation_x[i]*deviation_y[i] for i in range(len(x))]
    return sum(product_of_deviations)/len(x)
    
def correlation_coefficient(x,y):
    return bivariate_analysis(x,y)/deviation_value(x)/deviation_value(y)

####注意点

  • 標準偏差と共分散のプログラムを使っています。

##おまけ編
n人のクラスがあります。誕生日が同じの生徒が存在する確率を求めよ。

##<おまけ編の実装編>
同じ誕生日の二人組がいる確率についてをつかって実装していきたいと思います。同じ誕生日の二人組がいる確率についてより、求める確率は

1−_{365}P_n\div365^n

になる。
#####プログラム

from math import factorial
def permutations_count(n, r):#下の※のリンクから抜粋
    return factorial(n) // factorial(n - r)
def birthday(n):#n人のクラスで同じ誕生日の人がいる確率
    return 1-permutations_count(365,n)/365**n

※Pythonで階乗、順列・組み合わせを計算、生成からP(n,r)のコードを抜粋した。

このプログラムを使って初めて75%を超えるのはnがいくつ以上の時であるか、求めてみましょう。また、その時の確率を求めてみましょう。

####プログラム

people = 0
while True:
    if birthday(people) >= 0.75:
        print(people)
        print(birthday(people))
        break
    people += 1

peopleはクラスの人数です。birthday(100)が75%を超えると言うことがわかっていれば、for文でも良いかもしれません。23人ですでに50.7%のようです。

##参考文献

5
4
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
5
4

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?