2022-05-03

統計学入門　第3章　練習問題　解答まとめ

統計学

3.1
3.2
3.3
3.4 ブーストラップ
- i)
- ii)
- iii)

3.1

　pythonで散布図を記述する．なお，データ量が多いため先頭10個に省略して記述している．

r=∑(xi−y¯¯¯)(yi−y¯¯¯)∑(xi−x¯¯¯)2−−−−−−−−−√∑(yi−y¯¯¯)2−−−−−−−−−√

\begin{align}
r = \frac{\displaystyle \sum(x_{i} - \overline{y})(y_{i} - \overline{y})}{\sqrt{\sum(x_{i} - \overline{x})^{2}}\sqrt{\sum(y_{i} - \overline{y})^{2}}}
\end{align}

であるため，定義通り計算し

となる．

3.2

　A氏：喫煙者にはストレス解消の為に吸っている人が多い．タバコと肺がんに関係があるので無くストレスと肺がんに関係があるのではないか．私はストレスを抱えてないため肺がんにはなりにくいはずだ．

　コメント：A氏の主張はタバコと肺がんは見かけ上の相関であるという主張である．これを反論するにはタバコとストレスの関係，ストレスと肺がんの関係を調べたり，偏相関係数を求める必要がある．

3.3

　(2)の短期大学生と(3)の経営者団体で求めてみる．

　スピアマンの順位相関係数は

\begin{align}
r_{s} = 1 - \frac{6}{n^{3}-n} \sum (R_{i} - R_{i}')^{2}
\end{align}

で定義される．よって，

となる．ケンドールの順位相関係数は

\begin{align}
r_{k} = \frac{G - H}{n(n-1)/2}
\end{align}

である．(G, Hについては第三章まとめを参照)よって，

となる．

3.4 ブーストラップ

　標本集団から標本集団と同じ数だけランダムに値を再抽出し，そのデータセットに対し統計量を求め，それを何回も繰り返すことで母集団の性質を求めようとする手法をブートストラップと言う．

i)

　pythonで整数の乱数はrandomモジュールのrandint関数で発生できる

ii)

　今後のため相関係数を求める関数を用意しておく．

この関数を使い，11個のランダムな番号からデータを抽出し，その相関係数を求めると

となる．

iii)

　ii)を200回繰り返し，得られた相関係数200個をヒストグラムにする．コードは

である．得られたヒストグラムは

である．この統計手法をブーストラップという．

2022-05-02

機械学習を解釈する技術　4章　まとめ

機械学習

PD
実データでの分析

PD

　PFIは特徴量をシャッフルすることで，その特徴量がどれほど予測値に影響を与えているのかを調べることが出来た．しかし，PFIは特徴量の重要度しか分からず，特徴量を上げると予測値が上がるのかどうかなど，特徴量と予測値の関係までは分からない．

　特徴量と予測値の関係を知ることが出来るのがPDである．PDでは，他の特徴量を固定して興味のある特徴量だけを動かし，各データの予測値を平均する．例えば，ある特徴量のみを増加し全てのデータに対し予測させ予測値の平均を取り，グラフにプロットすればその特徴量と予測値の関係が可視化されることになる．

　PDは複雑なブラックボックスの沢山の特徴量と予測値の関係を直接解釈するのは難しいが，一つの特徴量と予測値の関係に絞れば人間にも解釈できるだろうというアイデアから生まれている．

　学習済みモデルを $\widehat f(x_0, x_1)$ とする．今特徴量 $x_0$ と予測値の関係を知りたい場合，数式で表すと

\begin{align}
\widehat{PD}_0(x_0) = \frac{1}{N}\sum \widehat f(x_0, x_{i, 1})
\end{align}

となる．この数式を観察すると，PDは $x_1$ で周辺化し平均をとったものと言える．

　PFIと同様に，PDもあくまで特徴量と予測値の関係を表しているだけで，特徴量と予測値の因果関係を表している訳では無いことに注意する．

実データでの分析

使用データの確認

　今回も3章で使ったデータと同じspaceship titanicコンペのデータを使う．

モデルの学習と特徴量の重要度計算

　3章と同じ方法で学習し，特徴量の重要度を求める．(シード値を固定していないため)3章と微妙に結果が異なり最も重要度が高いのはFoodCourtという特徴量であった．

　今回はFoodCourtという特徴量のPDを求めてみることにする．

PDの実装

　前回はスクラッチで実装したが，今回はsklearnにあるライブラリを利用してみる．sklearnのinspectionライブラリにはPDを計算するpartial_dependenceという関数がある．

scikit-learn.org

これを使うと今回の場合，

とたった2行で求めることが出来る．引数estimatorには学習済みモデル，Xには用いるデータ，featuresにはPDを求めたい特徴量，kindには予測値の平均を返してほしいならaverage，平均を取らないならindividualを入れる．(individualの場合ICEと呼ばれる指標になる．これは5章で解説される)今回はFoodCourtのPDを計算したいのでfeaturesにFoodCourtを入れた．

　この関数は予測値の平均と特徴量の値をdictで返す．これをグラフにプロットための関数もsklearnには存在する．

scikit-learn.org

　コードは以下になる．

　本書に従って新たに関数化しているが，肝心のコードは5行目の1行分のみである．取る引数はほとんど変わらない事が分かる．プロットされるグラフはこのようになった．

　X軸にある太い線は元データがどこに密集しているかを表している．FoodCourtは宇宙船内で請求された食事に関する金額である．よって元データでは請求された金額が約0，250，1100ドルの乗客が多いことが分かる．

　グラフを見るとFoodCourtが大きくなるほどPDが低くなっていることが分かる．つまり，生き残る確率が低くなっていると言える．請求された金額が高い＝沢山食事を取った乗客は生き残りにくそうだと言えるので，このモデルの結果に問題はなさそうだと考えられる．

2022-04-30

Atcoder 灰diff埋め ABC083C~100C 振り返り

競技プログラミング

灰diff埋めは簡単な問題が多いためテンポ良く大量に解いていこうと思っていましたが初期の問題は意外に難しい問題が多く……

ABC086 C - Traveling【AC】

　迷路の問題を解くように幅優先探索で解こうとすると，t=10^5もあるため間に合いません．この移動は同じ場所に留まる事が出来ない事を利用して，各時刻に対してO(1)で求めることが出来ます．
　時刻0で(0, 0)におり，時刻が1増えるごとにxかyのどちらかに-1か+1移動します．つまり，時刻1では(1, 0), (-1, 0), (0, 1), (0, -1)のどれかにいることになります．このことから，時刻の偶奇とx+yの偶奇が一致することが分かります．時刻1ではどう頑張ってもx=2の所に移動出来ませんし，時刻2ではどう頑張っても(1, 0)に移動出来ません（同じ場所に留まる事が出来ないため）
　そのため，まず時刻の偶奇とx+yの偶奇をチェックすれば良いです．しかし，例えば時刻2で(100, 100)という場合このチェックでは正しくありません．偶奇は合っていても距離が届かない場合をチェックする必要があります．これは現在の時刻とx, yの地点の差を求めれば良いです．x, y地点の方が現在の時刻より多ければ，最短で移動しても届かないのでNoとなります．

pre_t = 0
pre_x = 0
pre_y = 0
for _ in range(int(input())):
    t, x, y = map(int, input().split())
    if t % 2 != (x+y) % 2 or t - pre_t < abs(x - pre_x) + abs(y - pre_y) :
        print("No")
        break
    pre_t = t
    pre_x = x
    pre_y = y
else:
    print("Yes")

ABC100 C - *3 or /2【AC】

　最大で何回できるか求める必要があります．3倍すると言う操作は無限にする事が出来るので考える必要がありません．必ずどこかの数字は2で割る必要があり，数字は整数である必要があるため，なるべく2で割る操作を少なくすることが最大回数を求めることになります．なるべく少なくするので，1回の操作で2で割る数字は一つのみで他は3倍する操作をした方が良いです．結果，各数字で何回2で割る事が出来るかを数えればよいことが分かります．もし3倍で無く偶数倍だったら，無限に操作ができてしまいます．

N = int(input())
A = [int(i) for i in input().split()]
ans = 0
for a in A:
    while a % 2 == 0:
        a /= 2
        ans += 1
print(ans)

2022-04-28

Atcoder 灰diff埋め ABC079C 振り返り

競技プログラミング

atcoder problemsを見て，灰diffで未ACの問題を解いていきます．灰diffは簡単な問題も多いので，一部を抜粋して

ABC079 C - Train Ticket【AC】

　見るからに冗長なコードになっているのが残念です．

S = input()
A, B, C, D = int(S[0]), int(S[1]), int(S[2]), int(S[3])

for i in range(2**3):
    ans = A
    anslis = []
    if (i >> 0) & 1:
        ans += B
        anslis.append("+")
    else:
        ans -= B
        anslis.append("-")
    if (i >> 1) & 1:
        ans += C
        anslis.append("+")
    else:
        ans -= C
        anslis.append("-")
    if (i >> 2) & 1:
        ans += D
        anslis.append("+")
    else:
        ans -= D
        anslis.append("-")

    if ans == 7:
        print(f"{A}{anslis[0]}{B}{anslis[1]}{C}{anslis[2]}{D}=7")
        break

　AC後に全ての提出からこの問題を解いた他のコードを見てみました．パッと見た感じ一番シンプルに書かれていたのがこのコードです．

import itertools
def find_eq(abcd):
    a, b, c, d = abcd
    ops = ("+", "-")
    for op1, op2, op3 in itertools.product(ops, repeat=3):
        eq = a + op1 + b + op2 + c + op3 + d
        if eval(eq) == 7:
            return eq + "=7"

　このコードから学んだ事はいくつもあります．まず，

a,b,c,d = x

と書くことで文字列を自動的にアンパック(?)する事が出来るということです．シンプルにこう書くことで実装出来るということは初めて知りました．

x = ("+", "-")
for i, j, k in itertools.product(x, repeat=3):
    print(i, j, k)

>> +++
>> ++-
>> +-+
>> ...

と-の組み合わせをどう記述するかが一番の悩みどころでしたがitertoolsを使う事でこんなにシンプルに書けるとは知りませんでした．

eval("1+1")
>> 2

演算子+, -をどうやって式に組み込むか分からずにif文をダラダラ書く事になってしまいましたが，eval関数を使えば良いと知りました．eval関数は引数を式として評価するらしく，文字列を引数に入れても式として評価され数字の結果が返ってきます．
この問題はdiff337でしたが，昔の問題と言うこともあり学ぶ所が多く良かったです．

2022-04-27

Atcoder ABC249 振り返り

競技プログラミング

　最近土日に用事が重なってリアルタイム参加出来てません……

A - Jogging【AC】

　毎秒進むか停止するかを判断し，X秒までシミュレーションして進んだ距離が大きい方を表示します．

A, B, C, D, E, F, X = map(int, input().split())
takahashi = 0
aoki = 0

for i in range(X):
    if i % (A+C) < A:
        takahashi += B
    if i % (D+F) < D:
        aoki += E

if takahashi > aoki:
    print("Takahashi")
elif takahashi < aoki:
    print("Aoki")
else:
    print("Draw")

B - Perfect String【AC】

　全ての文字が相異なるという条件があります．これは，文字列をsetにとり元の文字列と長さを比較すれば良いです．同じ文字列があったらsetでは同じ文字列は消されるので元の文字列より短くなります．
　次に大文字と小文字の判定をします．upperとlowerというフラグ変数を用意し，始め両方Falseにしておきます．大文字があったらupper=True, 小文字があったらlower=Trueとします．大文字しか無かったり小文字しかない文字列だったらどちらかがFalseのままなのでそこで判定すれば良いです．

S = list(input())
if len(set(S)) != len(S):
    print("No")
else:
    upper = False
    lower = False
    for s in S:
        if s.isupper():
            upper = True
        if s.islower():
            lower = True

    if upper and lower:
        print("Yes")
    else:
        print("No")

C - Just K【AC】

　Nが15と小さいため全探索しても間に合うと判断出来ます．bit全列挙を使うとシンプルに書くことが出来ます．
　0か1がN個並んだ数字について，iが1ならi番目の文字列を選ぶ，iが0なら選ばないというようにします．このとき0から2^Nまでの数字を2進数化すると，001, 010, 011, 100, 101, 110, 111のように，選ぶ選ばないのパターンを全て列挙することが出来ます．
　i >> jというのは，iをjビット右に移動するという意味です．例えば，3 >> 1なら011を1ビット右に移動するので01となります．i & jはandです．iが1かつjが1の時に1，それ以外の時に0を返します．iが何行もある場合，一番右端のビットについてandが計算されます．
　つまり，(i >> j) & 1はiのj番目のビットが1なら1を，0なら0を返す処理になります．j番目が1なら文字列を選ぶので，リストに入れるというような処理を書けばいいです．
　ビットによって選ぶ文字列を選択出来たら，各文字が何個あるかを辞書で数えて，Kと等しいものの数を解とします．各文字列最大でも26文字しかないので，s.count(s_list)==Kでも判定できるかなと思いましたがTLEとなりました．

N, K = map(int, input().split())
S_list = [input() for _ in range(N)]
ans = 0

for i in range(2**N):
    s_list = []
    for j in range(N):
        if (i >> j) & 1 == 1:
            s_list += list(S_list[j])

    ans_cand = 0
    dic = {}
    for s in s_list:
        dic.setdefault(ord(s), 0)
        dic[ord(s)] += 1

    for v in dic.values():
        if v == K:
            ans_cand += 1

    ans = max(ans, ans_cand)

print(ans)

2022-04-26

機械学習を解釈する技術　3章　まとめ

機械学習

PFI
LOCOFI
GPFI
因果関係
PFIの実装と実データでの分析

PFI

　線形回帰なら回帰係数の大きさを見ることで特徴量の重要度を簡単に確認出来る．しかし，RandomForestなどの決定木モデルには回帰係数に相当する物が無く，NeuralNetはパラメータが複雑に関わり合い，パラメータから解釈するのは困難である．よって，どんなモデルでも計算できるアプローチを取る必要がある．

　どんなモデルでも共通して持っている要素は，予測値を出力する事である．よって，ひとつひとつの特徴量に対し，その特徴量が使えない時の予測値を計算して，全ての特徴量が使えるときの予測値と比較し，その差を見るという手法が考えられる．

　差が大きいと言うことは，その特徴量を使う事で性能が向上する＝重要であると言えるし，差が小さかったらその特徴量を使わなくても似た性能を出せる＝重要でないと言える．このように，予測値の誤差を調べることで特徴量の重要度を定義する手法がPFIである．

　PFIは特徴量の値をシャッフルすることで，特徴量を意味の無い状態＝使えないようにする手法である．ある特徴量ひとつをシャッフルさせたデータで予測値を出し，シャッフルさせてないデータで予測値を出し，両者の差を特徴量重要度とする．これを全ての特徴量に対して行うことで，特徴量の重要度を求めることが出来る．

LOCOFI

　PFIはデータをシャッフルさせることで特徴量を使えなくしていたが，LOCOFIではより直感的に，データを削除する事で特徴量を使えなくする．ただし，特徴量を削除するとモデルの構造も変わってしまうため，学習もし直さないといけない．さらに，モデルの構造が変わってしまうため，本来評価したいモデルと同じ特徴量重要度を求めることが出来ているとはいいがたいという欠点もある．

GPFI

　特徴量同士が強い相関関係にある場合，PFIだと重要度が低くでる事がある．強い相関を持つ二つの特徴量があったとき，ある一方の特徴量を使えなくしても，相関関係のあるもう一方の特徴量を使えばある程度高い性能の予測値を出すことが出来てしまうからである．　

　このような場合，相関関係の強い特徴量のグループはひとつにまとめてシャッフルしてしまい，一つの特徴量として重要度を出すという手法が考えられる．これがGPFIである．

　相関関係だけでなく，one hot encodingした特徴量でも重要度が正しく求められない事がある．one hot encodingとは，例えば出身地が日本，アメリカ，イギリスのどれかを取るデータがあったときに，日本という特徴量を新たに作り，出身地が日本なら1そうでないなら0とする手法である．同様にアメリカという特徴量を作り，アメリカ出身なら1そうでないなら0とする．

　このような特徴量についてPFIを行った場合を考える．例えば日本の重要度を求めるためにシャッフルする．すると，日本が0だったデータが1となる可能性がある．しかし，そのデータはアメリカかイギリスで1のはずで，日本でも1となると出身地が日本とアメリカというようにあり得ない事になる．このようなデータを使って予測すると正しい重要度が得られない可能性がある．

　このような問題を回避するために，GPFIでカテゴリカル変数をまとめてシャッフルしてしまうという方法がある．まとめてシャッフルする事で出身地が複数あるような矛盾を回避することが出来る．

　また他に，label encoding使うという方法がある．one hot encodingのように，日本アメリカイギリスという特徴量を作るのではなく，出身地という特徴量を作り，日本なら1アメリカなら2イギリスなら3と指定する方法である．これならシャッフルしても日本とアメリカが出身地になるような矛盾は生じず，正しい重要度が得られると思われる．

　一般に，one hot encodingの方が高い性能を出す事が多いが，カテゴリカル変数を使う事で重要度の解釈性が高まる．逆に，one hot encodingを使ってPFIやGPFIを使ったときは，注意する必要がある．

因果関係

　PFIで求まった特徴量の重要度を因果関係として解釈する事は良くないとされている．理由の一つに疑似相関の可能性があるからである．

　因果関係のあるA→Bと，相関関係のあるA, Cがあったとき，CとBには高い関係が見られる．しかしC→Bという因果関係は無いという事は十分ありうることであり，これを疑似相関と呼ぶ．

　PFIで求まった特徴量の重要度を見て疑似相関か否かを判断するのは不可能なので，重要度から因果を見いだしてはいけない．

PFIの実装と実データでの分析

　PFIを実装し，kaggleにあるspaceship titanicコンペのデータを使って動作を確認してみる．（タイタニックコンペのデータを使おうと思っていましたがその宇宙版のコンペが新しく追加されていたのでこちらにしてみました）

使用データの確認

Spaceship Titanic | Kaggle

からデータをダウンロードします．kaggle APIを使うと便利です．

　train.csvをpandasで読み込んでデータの一部を見てみます．

出身地やCabinの情報などいろいろなデータがある事が分かる．（全体的にタイタニックコンペと似ていますがHomeplanetなど宇宙チックでユニークなデータがあります）

モデルの学習

　このコンペは転送されたかされてないか(生き残ったか生き残ってないか？)を予測するコンペなので分類タスクになる．よって使うランダムフォレストも分類のものを使う事に注意する．ランダムフォレストは数値データしか使えないため，カテゴリカル変数をlabel encodingする．すると以下になる．(今回は特徴量の重要度を見るための実験なので特徴量エンジニアリングなどの工夫はなしでシンプルにいきます)

PFIの実装

　基本は本書の実装に沿います．

特徴量の重要度の表示

　特徴量の重要度を表示します．

　乗客が航海中に仮死状態になることを選択したかどうかを示す(!?)CryoSleepという変数や，食事に関して宇宙船で請求された金額を示すFoodCountが特徴量の上位に来ていることが分かります．

　PFIを用いることでSpaceship Titanicのデータの特徴量の重要度を求めることが出来た．

2022-04-25

統計学入門　まとめ

統計学まとめのまとめ

統計学入門 (基礎統計学Ⅰ)

東京大学出版会

Amazon

まとめ
各章まとめ
- 第3章
- 第4章
- 第5章
- 第6章
- 第7章
- 第8章
- 第9章
- 第10章
- 第11章
- 第12章

まとめ

　統計検定2級を取得しましたが，今一度統計の基礎をしっかり押えたいと思い本書の勉強を始めました．レビューなどでは数学的に難易度が高いという意見が多く，躊躇してましたが，実際に読んでみると大体の所は躓くことなく読み進めることが出来ました(入門というには難しいと思いますが……)．

　本書で出てくる各式に対して，自分がすぐ理解出来るレベルまで途中式や変形を略すことなく記述することが出来ました．統計検定準1級の取得を目標に頑張っていきます．

　次は本書の練習問題を解いていきます！

各章まとめ

3.1

3.2

3.3

3.4 ブーストラップ

i)

ii)

iii)

PD

実データでの分析

使用データの確認

モデルの学習と特徴量の重要度計算

PDの実装

ABC086 C - Traveling【AC】

ABC100 C - *3 or /2【AC】

ABC079 C - Train Ticket【AC】

A - Jogging【AC】

B - Perfect String【AC】

C - Just K【AC】

LOCOFI

GPFI

因果関係

PFIの実装と実データでの分析

使用データの確認

モデルの学習

PFIの実装

特徴量の重要度の表示

まとめ

各章まとめ

第3章

第4章

第5章

第6章

第7章

第8章

第9章

第10章

第11章

第12章