乙Py先生のプログラミング教室
初学者のためのプログラミング学習サイト

データサイエンス

さて、今回は機械学習のお話です。

一応なんちゃって
MLエンジニアなのでねwwwwwww

食べログさんのお店をデータ化したら
星の予測モデルが作れました。

解説動画はこちら





以前の動画では
食べログ分布というものについて触れました。
前回の記事はこちら食べログ分布は有るのか

今回は
実際に星を予測するモデルを作っていきましょう。

さて
こういった機械学習を用いた
予測モデルを作成するには
まずはデータが必要です。


こちらは
みなさん頑張って集めてみて下さい!!!!!

データとしては
こんなデータです。

スクリーンショット 2019-10-26 21.13.18

コメント数、ブックマーク数や価格
コメントを残した人の星の平均などです。

そして予測をするのは
「星」のデータとなります。


予測モデルの作成に使用したコードは
次のようなコードになります。
import warnings
warnings.filterwarnings('ignore')

import pandas as pd
import numpy as np

from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
import xgboost as xgb
import lightgbm as lgm

from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_squared_error

import matplotlib.pyplot as plt
%matplotlib inline


data_df = pd.read_csv('tabelog_star_data.csv')

data_df = data_df[data_df['星']>2.9]

X = data_df.drop(['星'],axis=1)
Y = data_df['星']
x_train, x_test, y_train, y_test = train_test_split(X,Y, test_size=0.3)

models = []
models.append(("DesicionTree",DecisionTreeRegressor(max_depth=8)))
models.append(("RandomForest",RandomForestRegressor(max_depth=8)))
models.append(("LightGBM",lgm.LGBMRegressor(max_depth=8)))
models.append(("XGBoost",xgb.XGBRegressor(max_depth=8,objective='reg:squarederror')))

result = {}
for name,model in models:
    print(name)
    model.fit(x_train, y_train)
    y_predict = model.predict(x_test)
    mae  = mean_absolute_error(y_test, y_predict)
    rmse = np.sqrt(mean_squared_error(y_test, y_predict))
    result[name] = [mae,rmse]
    print('  MAE : ',mae)
    print('RMSE : ',rmse)
    print()

少しづつみていきましょう。

まずは
データを読み込みします。

pd.read_csv()でCSVファイルを読み込みして
データフレームに変換できます。

なおこのデータは
すでに「前処理」と呼ばれる
機械学習用のデータに
整形する工程を終えたデータで
全て数値型のデータになっています。


次に訓練用とテスト用のデータに分けます。

train_test_splitを用いると
簡単に振り分けることができるので楽です。

今回は3割テスト、残りを訓練用としました。

モデルは
今回は自分がよく使っている
決定木,ランダムフォレスト,LightGBM,XGBoost

こちらの4種類を使っていきます。

モデルの作成は
model.fit(x_train, y_train)

これだけで済んでしまいます。
機械学習と言っても
コードはたったこれだけなんですよね・・・

モデルの作成が終わったら
次はテストデータで予測を行い
その精度を検証します。

y_predict = model.predict(x_test)

予測もたったこれ1行だけで済んでしまいます。

あとはこの結果から
どれだけ誤差が出たのかを見て
精度を図ります。

誤差を見るのに
MAE
RMSE
を使います。

今回はざっくり検証なので
どれだけ誤差が出たのかだけを見ることとします。

DesicionTree
  MAE :  0.1047453451900971
RMSE :  0.14276028664764914

RandomForest
  MAE :  0.08806766700129451
RMSE :  0.11329995353232411

LightGBM
  MAE :  0.08365354791749746
RMSE :  0.10874470218772024

XGBoost
  MAE :  0.08733976372263648
RMSE :  0.11181978441280738

MAEはどんだけ誤差が出たかを平均したものです。
MAE0.083なら
+-0.083くらいはズレます。

実際の星が3.5だとして
大体予測が3.4-3.6くらいになるだろうという
モデルができました。

今回一番良かったものは
LightGBMですね。

誤差を可視化してみましょう。
tabe1

平均で0.083ズレですが
大きいところでは0.36くらいはズレてますね

誤差の累積だと
tabe2

80%くらいで0.15以内には収まりそうな感じですね。

予測の精度としては
まあまあなんじゃないでしょうか

これが仕事なら
クロスバリデーションしたり
パラメータいじったり
元々の特徴量作る集計の見直しなど
まだまだやること有ります。

色々すればもう少し
精度は上がるでしょうが
ざっくり検証なのでこれまでです。


さて
このデータのうちで
どれが一番予測に寄与しているでしょうか?

寄与度を出してみました。


191コメント数

181ブックマーク数

1711rate_1000_s_mean

1521rate_100_s_mean

1491rate_5000_s_mean

1341rate_500_s_mean

106夜価格低

1051year_500_s_mean

811rate_2000_s_mean

721year_1000_s_mean

683month_500_s_mean

656month_500_s_mean

631rate_10000_s_mean

55昼価格低

こんな感じでした。

コメント数やブックマークの数
コメントした人の星の平均点
ディナー価格などが
関係しているようです。

実際の星の計算は
もう少し複雑な計算方法で
違うものでは有ると思いますが
近い形の星を算出できたと思います、


これで、ブラックボックスであった
星の計算方法に少し
近づいたんじゃないかなーと思います。

食べログに掲載されているお店さんで
点数を上げたければ
ここら辺がポイントになるのではと
思ってはおりますが・・・

このデータには
お店側が食べログ側に
お金を支払っているかどうか
というデータが含まれておりません。

なので、実はお金を支払っていないから
星が下がったということも
考えられなくは無いので
実際の予測にはこの情報が不可欠です。

どなたか
このフラグ値お持ちでは
無いでしょうかwwwwwwww

今回はここまでです
それでは。


先日
ブラックホールの撮影に成功したという
ニュースがありましたね。

なので
そのブラックホールの大きさを
プログラムで求めてみました。


解説動画はこちら



まずはブラックホールって
何なのって話ですが

ウィキペディアによると
極めて高密度かつ大質量で
強い重力のために物質だけでなく
光さえ脱出することができない天体である。

まあ元は星だったりとか
そんなやーつですね。

そのあとの説明では

周囲は非常に強い重力によって時空が著しくゆがめられ
ある半径より内側では脱出速度が光速を超えてしまう。

この半径をシュヴァルツシルト半径
この半径を持つ球面を事象の地平面(シュヴァルツシルト面)と呼ぶ。

つまりこの事象の地平面が
ブラックホールの大きさになるわけですが

その大きさを調べるには
シュヴァルツシルト半径が分かれば
その直径も分かりますね。

シュヴァルツシルト半径の
wikiの記事も見ると

1916年カール・シュヴァルツシルトは
アインシュタインの重力場方程式の解を求め
非常に小さく重い星があったとすると
その星の中心からのある半径の球面内では
光も脱出できなくなるほど曲がった時空領域が出現することに気づいた。

その半径をシュワルツシルト半径と呼び
シュワルツシルト半径よりも小さいサイズに収縮した天体は
ブラックホールと呼ばれる。

はい。

んでそれを求める公式があるそうです。

スクリーンショット 2019-04-16 21.51.54

これが関係する公式です。

早速いろいろ計算していきましょう。

まずは地球を縮めたら
どれくらいの大きさでブラックホールになるのかです。

重力定数と光速度は
scipyのconstants
というライブラリの中に定数が有るので
それを用いると計算がスムーズです。

まずライブラリの呼び出しは

from scipy import constants

お次は半径を求める式を
プログラムで当てはめていきます。

E = 5.9724 * (10**24)
G = constants.G
C = constants.c
R = 2 * G * E / (C**2)
print(R * 100)
0.8870107529844293

さて結果は・・・
0.88センチ

地球を1センチ程度に圧縮したら
ブラックホールになり得るのだそう。

次は太陽です。
S = 1.989 * (10**30)
G = constants.G
C = constants.c
R = 2 * G * S / (C**2)
print(R)
2954.0291803731

太陽は約3キロメートルほどで
ブラックホールになるんだそう。

逆に
半径から質量を求めてみましょう。
公式ではRを求めていたので
質量Mを求めるようにするには

R = 0.008870107529844293
M = R * (C**2) / (2 * G)
print(M)
5.9724e+24

はい
これで半径から質量も
出すことができましたね。

最後にこの前撮影された
ブラックホールの半径を求めます。

撮影できたものの質量は
太陽の65億倍なんだとか

この数値を使って当てはめてみると・・・

S = 1.989 * (10**30)
B = S * 6500000000
G = constants.G
C = constants.c
R = 2 * G * B / (C**2)
print(R/1000)
19201189672.42515


約192億キロです。

デカすぎワロタwwwww

よくそんなの撮れましたね。

引くほどデカイ!!!


ブラックホールといえば
あの名作
インターステラーではないでしょうか



滅亡しかけている地球のために
恒星間移動をして
住める惑星を探して調査して
一難去ってまた一難

一悶着ありまくりの後
最後の賭けで
ブラックホールに入って
この世界の謎を解き明かして
ラストは・・・

うーーん
奥が深かった。

理論物理学者のキップ・ソーンが
科学コンサルタント兼製作総指揮を務め
理論武装も完璧であろうと思われる。

物理学の話でもあるけど
親子の絆を描く作品でもある
不朽の名作

ぜひお試しあれ。


さて
またまた東大の問題です。

確率を求めたりする計算は
プログラムの得意とするところで
業務でもよく出てきます。

いかに東大の問題と言えども
プログラミングができれば
簡単に解けてしまうと思います。

解説動画はこちら




2019年の問題

正八角形の頂点を反時計回りにA,B,C,D,E,F,G,Hとする。
また投げたときにおもて裏の出る確率がそれぞれ2分の1のコインがある。

点Pが最初に点Aにある。次の操作を10回繰り返す。
操作:コインを投げ表が出れば点Pを反時計回りに隣接する頂点に移動させ
裏が出れば点Pを時計回りに隣接する頂点に移動させる。
例えば点Pが点Hにある状態で、投げたコインの表が出れば点Aに移動させ
裏が出れば点Gに移動させる。

以下の事象を考える。

事象S:操作を10回行った後に点Pが点A上にある。
事象T:1回目から10回目の操作によって、点Pは少なくとも1回点Fに移動する。
(1)事象Sが起こる確率を求めよ
(2)事象Sと事象Tがともに起こる確率を求めよ


まずは問題文だけでは分かりづらいので
可視化してみましょう。

正8角形
反時計回りにA,B,C,D,E,F,G,H

この二つの条件で図形を描画してみます。

ライブラリを読み込んで
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
fig = plt.figure(figsize=(2,2),dpi=200)

r = 8
x = [np.sin(i) for i in np.linspace(0,2*np.pi,r+1)]
y = [np.cos(i) for i in np.linspace(0,2*np.pi,r+1)]
plt.plot(x,y)

t = ['A','H','G','F','E','D','C','B','']
for i ,(x1,y1) in enumerate(zip(x,y)):
    plt.text(x1,y1,t[i],ha='center',va='bottom')

plt.show()
8kakkei

このように描画できます。

x,yには正8角形の各x,yの座標

plt.plotではx,y座標を指定し
plt.textで図形に文字を付け加えられます。

生成した座標値は時計回りの座標なので
それを見越して表示させる文字を変えています。

それではイメージできたところで
問題を解いていきましょう。

1つめの問いは

(1)事象Sが起こる確率を求めよ
事象S:操作を10回行った後に点Pが点A上にある。 

まずは問題を解くのに必要なライブラリを読み込みます。
import itertools
from fractions import Fraction
coins = [i for i in itertools.product([1,-1],repeat=10)]
print(len(coins))

base = ['C','B','A','H','G','F','E','D','C','B','A','H','G','F','E','D','C','B','A','H','G’]
print(base[10])
1024
A

コインを10回ふった際の結果をcoinsという変数に格納します。

itertools.productで表裏の結果10回繰り返した際の順列が生成できます。
引数はその1回の組み合わせをリストで定義し
repeat=回数 で繰り返し回数が指定できます。

表は-1、裏は1として定義してrepaetが10回ですね。

単純に2の10乗回数だけコインをふることになりますね。

baseという変数にはAから始まって
時計回りに最大10回
反時計回りに最大10回
移動した際の点を入れます。
真ん中のbase[10]が開始地点です。


用意ができたらあとは計算だけですね。
10回ふった際の最終地点を求めて
その結果を辞書に格納してあげます。

最終地点:回数

とすることで確率が求められます。

calc_dict = {}
for coin in coins:
    key = base[10+sum(coin)]
    if key in calc_dict:
        calc_dict[key] +=1
    else:
        calc_dict[key] =1

for k,v in calc_dict.items():
    print(k,v)
C 256
A 272
E 240
G 256

print(calc_dict['A']/len(coins))
print(Fraction(calc_dict['A'],len(coins)))
0.265625
17/64

コインを10回ふった組み合わせ1024回のうち
Aに戻るのは272回でした。

なので確率としては
272/1064 

約分して
17/64となります。

Fractionで分数計算ができて
約分もしてくれます。

Pythonはこういうところが
めちゃくちゃ便利です。


2つめの問いは

(2)事象Sと事象Tがともに起こる確率を求めよ
事象S:操作を10回行った後に点Pが点A上にある。 
事象T:1回目から10回目の操作によって、点Pは少なくとも1回点Fに移動する。

条件が二つに増えます。

それではコインを10回ふった際の
軌跡を求めてみましょう。
coin_result = []
for coin in coins:
    tmp = 10
    result = ''
    for c in coin:
        tmp += c
        result += base[tmp]
    coin_result.append(result)

組み合わせの数だけ試行して
移動結果を文字としてつなげます。
リスト型に結果を保存します。

最終的に求める条件は
1.Fを通る(Fがある)
2.Aで終わる
この二つを満たす結果があれば
1カウントし、全組み合わせの中から
合致するものの個数を数えます。
count = 0
for row in coin_result:
    if 'F' in row and row[-1]=='A':
        count+=1
print(count)

66
print(count/len(coins))
print(Fraction(count,len(coins)))
0.064453125
33/512

最終的な答えの確率は
66/1024

約分して
33/512
となりました。

さていかがだったでしょうか?
確率を求めたりする計算は
日常業務ではよくでてきます。

数学的な可視化や
その計算については
numpyやmatplotlibライブラリを
用いることが多いです。

単純な計算については
通常のPythonプログラムだけでも
行うことができます。

問題を読み間違わなければ
このくらいの計算は
プログラミングで簡単に行うことができるので

プログラミングを覚えていない方
これから覚えたい方は
ぜひPythonを覚えてみてください。

Pythonについては無料の動画講座を用意しています。
乙py式5時間で学ぶプログラミング基礎(python編)

興味のある方はぜひこちらをご参照くださいませ。

それでは。


今回はヒマなので
なーんか面白いことないかーとか
思っていましたら

東の方に
有名な大学が有るらしいんですよ。


東京大学って知ってます?

自分は行ったことないんで
よく知らないんですけど
有名大学みたいなので
数学の問題を解いてみることにしました。

解説動画はこちら



2003年の東大の入試問題より
 
円周率が 3.05 より大きいことを証明せよ。


さて
どう解きましょうかねー

とりあえず
円周の長さ > 内接X角形の外周

になるはずなので
この内接する多角形の外周を求めて
それが3.05よりも大きければ

総じて
円周率> 3.05 
となるのではないかと思います。

まずは円を描いてみましょう。

Pythonのライブラリを用いて
円を描いてみます。

# ライブラリの読み込み
import numpy as np
import matplotlib.pyplot as plt
% matplotlib inline


描画用のライブラリとして
matplotlibを読み込んでおきます。

# 半径5の円を描く
plt.figure(figsize=(2,2),dpi=300)
r = 5
x = [np.sin(np.radians(_x))*r for _x in np.linspace(-180,180,361)]
y = [np.cos(np.radians(_y))*r for _y in np.linspace(-180,180,361)]
plt.plot(x,y)

plt.xticks([i for i in range(-r,r+1)])
plt.yticks([i for i in range(-r,r+1)])
plt.axes().set_aspect('equal','datalim')
plt.show()
en1
円を描くにはx,y座標が必要ですが
その点を求めるのにnumpyを使います。

numpy.linspaceで等間隔の配列が作成できます。

あとはx,y座標の点を求めるのに
numpy.sin , numpy.cos を使います。
sin , cos に渡せるのはradianでないといけないので
一旦numpy.radiansで変換します。

これで円を描くことができます。


次にこの円に内接する多角形を描いてみましょう。
今回は12角形を描くこととします。

わかりやすくするために
三角形の辺の長さが3:4:5になるという法則を用いて
座標を決めていきます。

一番長い5の辺は円に接する先端ですね。
あとはx,yの座標は3か4になります。

一番上をx=0,y=5として座標点をプロットしてあげます。
# 内接する12角形を描画
plt.figure(figsize=(2,2),dpi=300)
r = 5
x = [np.sin(np.radians(_x))*r for _x in np.linspace(-180,180,361)]
y = [np.cos(np.radians(_y))*r for _y in np.linspace(-180,180,361)]
plt.plot(x,y)

x2 = [0,3,4,5,4,3,0,-3,-4,-5,-4,-3,0]
y2 = [5,4,3,0,-3,-4,-5,-4,-3,0,3,4,5]
plt.plot(x2,y2)

plt.xticks([i for i in range(-r,r+1)])
plt.yticks([i for i in range(-r,r+1)])
plt.axes().set_aspect('equal','datalim')
plt.show()
en2
はいこれで円に接する12角形が描けました。

あとはこの辺の長さを求めて
あげれば良いということになります。

辺の長さを求めるにはどうすれば良いでしょうか?

これはnumpyを使って
ユークリッド距離を求めることで
辺の長さを計算することができます。

numpy.linalg.norm(座標a - 座標b)

これで2点のユークリッド距離を求めることができます。

12角形のうち
全部を求める必要はなく

右上部分の4点を用いて
3辺の長さを計算してみましょう。

x,y座標はそれぞれ
0,5
3,4
4,3
5,0
となるので
そのうち2点を使って計算します。
l = np.array([[0,5],[3,4],[4,3],[5,0]])
c1 = np.linalg.norm(l[0]-l[1])
c2 = np.linalg.norm(l[1]-l[2])
c3 = np.linalg.norm(l[2]-l[3])
ans1 = 4 *(c1+c2+c3)
print(ans1)

30.9550755308

さてこれで多角形の辺の長さが計算できました。

問題文は3.05より大きいことを証明せよなので
比率を合わせます。

この円は半径5の円なので直径は10です。
なので10倍します。
30.5


10倍した30.5よりも
12角形の外周は30.95のため大きくなり

必然的にそれよりも円周は大きいので

円周率>3.05

になるはずです。

さて
描画などについては
Pythonのライブラリを用いると簡単に
描くことができ
また座標間の距離なども
簡単に計算することができます。

数学の問題では
このnumpyとmatplotlibライブラリを使って
いろいろ問題に応用することができます。

数学的な可視化や
その計算については
numpyやmatplotlibライブラリを
用いることが多いです。

単純な計算については
通常のPythonプログラムだけでも
行うことができます。

問題を読み間違わなければ
このくらいの計算は
プログラミングで簡単に行うことができるので

プログラミングを覚えていない方
これから覚えたい方は
ぜひPythonを覚えてみてください。

Pythonについては無料の動画講座を用意しています。
乙py式5時間で学ぶプログラミング基礎(python編)

興味のある方はぜひこちらをご参照くださいませ。

それでは。

さて
今回は確率のお話です。


「ロト7」の1等、同じ売り場から3口も当選

なんてニュースが最近有ったので
いろいろなギャンブルの確率が
どのくらいなのかを求めてみたくなりました。

解説動画はこちら



まずは競馬からです。

競馬は16-18頭の馬の中から
着順を当てる仕組みです。

馬券は正式には
勝ち馬投票券と言うらしいです。

16頭18頭だてでそれぞれを求めました。

確率計算にはPythonを用いますので
ライブラリをインポートします。

import math
import itertools
こちらを用います。

さて、まずは単勝から
単に1着を当てるもの

確率は
# 16頭
print('{0}/{1}'.format(1, 16))
print('{0:.5}%'.format(1/16*100))
1/16
6.25%
# 18頭
print('{0}/{1}'.format(1, 18))
print('{0:.5}%'.format(1/18*100))
1/18
5.5556%

次は複勝
1-3着に入る馬を当てます。

予想した馬が1位でも3位でも良いので
単純に単勝の3倍になります。
# 16頭
print('{0}/{1}'.format(3, 16))
print('{0:.5}%'.format(3/16*100))
3/16
18.75%
# 18頭
print('{0}/{1}'.format(3, 18))
print('{0:.5}%'.format(3/18*100))
3/18
16.667%


複勝であれば
6回に1回は当たるという確率ですね。

ただし払い戻しは平均2倍程度なので
あまり美味しくはありません。

お次は枠連
18頭の馬は9の枠に入れられます。
それぞれ1-2頭ほどで構成され
その枠での順位を当てるというものです。

9枠の中から2枠を選ぶというものです。

Pythonでの
組み合わせの求め方は2通りありますが
len(list(itertools.combinations(range(9),2)))

math.factorial(9) // (math.factorial(9 - 2) * math.factorial(2))
どちらであっても36通りです。

ということで確率は
waku = math.factorial(9) // (math.factorial(9 - 2) * math.factorial(2))
print('{0}/{1}'.format(1, waku))
print('{0:.5}%'.format(1/waku*100))
1/36
2.7778%

となります。


お次は馬単と馬連
選んだ2頭が1-2着になる組み合わせで
馬単は1-2着をそのまま当てる
馬連は1-2,2-1どちらでも良いという買い方です。

当然
馬連は馬単の2倍の確率になります。

馬練:
# 16頭
umaren= math.factorial(16) // (math.factorial(16 - 2) * math.factorial(2))
print('{0}/{1}'.format(1, umaren))
print('{0:.5}%'.format(1/umaren*100))
1/120
0.83333%
# 18頭
umaren= math.factorial(18) // (math.factorial(18 - 2) * math.factorial(2))
print('{0}/{1}'.format(1, umaren))
print('{0:.5}%'.format(1/umaren*100))
1/153
0.65359%

馬単:
# 16頭
umatan= math.factorial(16) // math.factorial(16 - 2)
print('{0}/{1}'.format(1, umatan))
print('{0:.5}%'.format(1/umatan*100))
1/240
0.41667%
# 18頭
umatan= math.factorial(18) // math.factorial(18 - 2)
print('{0}/{1}'.format(1, umatan))
print('{0:.5}%'.format(1/umatan*100))
1/306
0.3268%


お次はワイド
3着以内に入る組み合わせのうち2つを
それらの着順に依らず順不同で予想するものです。

1-2 , 1-3 , 2-3着になれば良いので
単純に馬連の3倍の確率になります。

# 16頭
wide = math.factorial(16) // (math.factorial(16 - 2) * math.factorial(2))
print('{0}/{1}'.format(3, wide))
print('{0:.5}%'.format(3/wide*100))
3/120
2.5%
# 18頭
wide = math.factorial(18) // (math.factorial(18 - 2) * math.factorial(2))
print('{0}/{1}'.format(3, wide))
print('{0:.5}%'.format(3/wide*100))
3/153
1.9608%


次は3連
3連単と3連複でそれぞれ
3連単:1着・2着・3着になる組み合わせを着順通りに予想する
3連複:1着・2着・3着になる組み合わせを順不同で予想する

3連複:は組み合わせ(16C3 , 18C3)
3連単:は順列(16P3 , 18P3)でもとめることになります。


3連複:
# 16頭
sanrenhuku = math.factorial(16) // (math.factorial(16 - 3) * math.factorial(3))
print('{0}/{1}'.format(1, sanrenhuku))
print('{0:.5}%'.format(1/sanrenhuku*100))
1/560
0.17857%
# 18頭
sanrenhuku = math.factorial(18) // (math.factorial(18 - 3) * math.factorial(3))
print('{0}/{1}'.format(1, sanrenhuku))
print('{0:.5}%'.format(1/sanrenhuku*100))
1/816
0.12255%


3連単:
# 16頭
sanrentan= math.factorial(16) // math.factorial(16 - 3)
print('{0}/{1}'.format(1, sanrentan))
print('{0:.5}%'.format(1/sanrentan*100))
1/3360
0.029762%
# 18頭
sanrentan= math.factorial(18) // math.factorial(18 - 3)
print('{0}/{1}'.format(1, sanrentan))
print('{0:.5}%'.format(1/sanrentan*100))
1/4896
0.020425%

このくらいの確率になると
でたらめに選んで買っても
5千回に1回しか当たりません。

そうとう当たらないので
当てるには絞り込みの技術が必要になってきます。

競馬の最後はWIN5
JRAが指定する5つのレース
それぞれの1着を予想するものです。

単勝の5乗ですね。

その確率は
win5 = 18**5
print('{0}/{1}'.format(1, win5))
print('{0:.10F}%'.format(1/win5*100))
1/1889568
0.0000529221%

約190万分の1の確率です。
当てた際の
払戻金の想定額が2億円ほどで

宝くじなどの確率が
1000万分の1の確率で
同じくらいの金額なので

WIN5の方が期待値は上かと思います。


競馬の確率の高さは

複勝>単勝>枠連>ワイド>馬連>馬単>3連複>3連単>WIN5

になります。

次は宝くじのうち
ロト6です。

ロト6は1-43まで数字のうち
本数字6個と1個のボーナス数字
計7個の数字を選択して当てるものです。

1-5等まであります。

まずは1等

本数字6個すべて一致

なので単純に
43C6 となります。
この数はロト6のベースとなる数です。

求め方は
loto_base = math.factorial(43) // (math.factorial(43 - 6) * math.factorial(6))
print('{0}/{1}'.format(1, loto_base))
print('{0:.10f}%'.format(1/loto_base*100))
1/6096454
0.0000164030%

1等は約600万分の1の確率となります。

2等以下は次のようになります。

2等 : 本数字5個と一致,ボーナス数字1個と一致
loto_base = math.factorial(43) // (math.factorial(43 - 6) * math.factorial(6))
print('{0}/{1}'.format(6, loto_base))
print('{0:.10f}%'.format(6/loto_base*100))
6/6096454
0.0000984179%

3等 : 6個のうち5個が本数字に一致
loto_base = math.factorial(43) // (math.factorial(43 - 6) * math.factorial(6))
loto3 = (math.factorial(6) // (math.factorial(6 - 5) * math.factorial(5)))*37-6
print('{0}/{1}'.format(loto3, loto_base))
print('{0:.10f}%'.format(loto3/loto_base*100))
216/6096454
0.0035430432%

4等 : 6個のうち4個が本数字に一致
loto_base = math.factorial(43) // (math.factorial(43 - 6) * math.factorial(6))
loto4 = (math.factorial(6) // (math.factorial(6 - 4) * math.factorial(4)))* \
((math.factorial(37) // (math.factorial(37 - 2) * math.factorial(2))))
print('{0}/{1}'.format(loto4, loto_base))
print('{0:.10f}%'.format(loto4/loto_base*100))
9990/6096454
0.1638657488%

5等 : 6個のうち3個が本数字に一致
loto_base = math.factorial(43) // (math.factorial(43 - 6) * math.factorial(6))
loto5 = (math.factorial(6) // (math.factorial(6 - 3) * math.factorial(3)))* \
((math.factorial(37) // (math.factorial(37 - 3) * math.factorial(3))))
print('{0}/{1}'.format(loto5, loto_base))
print('{0:.10f}%'.format(loto5/loto_base*100))
155400/6096454
2.5490227598%

さて
5等は1000円です。
2.5%の確率でしか当たらないので
40回買っても1000円しか戻ってこないのです。

1位がでなければキャリーオーバーで
持ち越しになり、金額が増えますが
そもそも600万分の1の確率では
人生において当たることは
ほぼないでしょうね。


今回のラストはポーカーです。

ポーカーはトランプ52枚を使ったゲームで
役が決まっています。

その役がどれだけ有るのかを求めます。

役は
一番強い役がロイヤルフラッシュ
(ロイヤルストレートフラッシュ)

ストレートフラッシュ: 連続した数字で(絵柄)
同じスートの5枚のカード。

ストレートフラッシュ同士の場合は、
最もランクが高いカードを持つプレイヤーが勝者
最も強いストレートフラッシュは同じスートのA-K-Q-J-10

フォーカード: 同じランクのカード4枚と他1枚カード
フォーカード同士の場合は、同じランクのカード
4 枚のランクが最も高いプレイヤーが勝者

フルハウス: 同じランクのカード 3 枚と別の同じランクのカード 2 枚。
(スリーカードとワンペアの組み合わせ)
フルハウス同士の場合は、3 枚 1 組のカードのランクが最も高いプレイヤーが勝者

フラッシュ: 同じスートの 5 枚のカード。
フラッシュ同士の場合は、最もランクが高いカードを持つプレイヤーが勝者

ストレート: 連続した数字の 5 枚のカード。
ストレート同士の場合は、最もランクが高いカードを持つプレイヤーが勝者
最も強いストレートは A-K-Q-J-T (エースハイ)、最も弱いストレートは 5-4-3-2-A (ファイブハイ)

スリーカード: 同じランクのカード 3 枚と、ランクの異なる 2 枚のサイドカード。
スリーカード同士の場合は、3 枚 1 組のカードのランクが最も高いプレイヤーが勝者

ツーペア: 同じランクのカード 2 枚 1 組が 2 組と、1 枚のサイドカード。
ツーペア同士の場合は、最もランクが高いペアを持つプレイヤーが勝者

ワンペア: 同じランクのカード 2 枚と、ランクの異なる 3 枚のサイドカード。
ワンペア同士の場合は、最もランクが高いペアを持つプレイヤーが勝者

ハイカード: 上記のどれにも当てはまらない手札。
ブタとも言います。

今回は役を判定するプログラムを作って
全カードの組み合わせからその役がなんなのかを判定して
役の組み合わせ数を求めます。

デッキはカード52枚として
それを判定する関数を作成し
全52枚の中から5枚を選んだ際の組み合わせから
役が何かを求めて役の個数をカウントします。
# デッキの生成
deck=[b+':'+str(a) for a in range(1,14) for b in ['C','S','D','H']]
# 役の判定
def jadge_role(card):
    s = {}
    for c in card:
        k = int(c.split(':')[1])
        if k in s:
            s[k]+=1
        else:
            s[k] =1
    t = {c.split(':')[0] for c in card}
    n = sorted([c for c in s.keys()])
    if len(t)==1 and all([1 in n,10 in n ,11 in n,12 in n,13 in n]):
        return 'RSF'
    if len(t)==1 and (all([1 in n,10 in n ,11 in n,12 in n,13 in n]) or
             (max(n)-min(n)==4) and len(s)==5):
        return 'SF'
    if 4 in s.values():
        return '4C'
    if 3 in s.values() and 2 in s.values():
        return 'FH'
    if len(t)==1:
        return 'FL'
    if (all([1 in n,10 in n ,11 in n,12 in n,13 in n]) or
             (max(n)-min(n)==4) and len(s)==5):
        return 'ST'
    if 3 in s.values():
        return '3C'
    if list(s.values()).count(2)==2:
        return '2P'
    if list(s.values()).count(2)==1:
        return '1P'
    return 'BT'
# 組み合わせを計算
cards = itertools.combinations(deck,5)
calc_dict = {}
for card in cards:
    role = jadge_role(card)
    if role in calc_dict:
        calc_dict[role] += 1
    else:
        calc_dict[role]   = 1

# 結果
poker_base = math.factorial(52) // (math.factorial(52 - 5) * math.factorial(5))
print(poker_base)
for k,v in sorted(calc_dict.items(),reverse=True,key=lambda x:x[1]):
    print(k,'\t',v,'\t','{0:.10f}%'.format(v/poker_base*100))
2598960

BT  1302540 50.1177394035%
1P  1098240 42.2569027611%
2P    123552 4.7539015606%
3C      54912 2.1128451381%
ST      10200 0.3924646782%
FL         5108 0.1965401545%
FH        3744 0.1440576230%
4C          624 0.0240096038%
SF            36 0.0013851695%
RSF           4 0.0001539077%

このような結果になります。
ロイヤルストレートフラッシュは
約64万回に1回ほどは出る計算ですね

フラッシュとストレートでは
ストレートの方が2倍も出る確率が高いのです。
なのでどちらかを迷う場合は
確率の高いストレートを狙うというのが
戦略になるのかもしれませんね。


いろいろなギャンブルにまつわる
組み合わせと確率を求めてみました。

次は
戦略なども考えていきたいかなと思っております。

それでは

このページのトップヘ