今回はゆるめで、あるイベントのLTに参加してきた感想やら、資料やら動画やらソースコードを共有するだけの記事となっております。

ことの成り行き

Tokyo.Rでお世話になっているy__mattuさんから、「Music×Analytics Meetup Vol.4」のイベントのことを教えていただき、せっかくなのでLTをさせてもらえることになりました。

音楽の分析に関して入門レベルなので果たして参加者が盛り上がってくれる何かを提供できるか不安でした。

納期的にも、今の実力的にも過去に書いたブログ([Python]スペクトグラムで楽器の試奏データを比較（プロ奏者とも）)の延長線上のことでいいなと思ったんですが、
先行研究は欲しいし、分析対象の音楽もスケール(音階)ではなく何かの曲が良いと思いました。あと、意外な情報も載せたい、分析結果を活かして便益を得たいとも思いました。

なので、LTまでに先行研究を見つけ、意外な情報としてjupyter notebookで録音する機能を見つけ、分析対象の曲にギャツビーのCMのものを採用しました。
そして、この研究から自分の楽器の演奏の腕前を向上させることができました。
これならLTとして発表しても入門レベルでも楽しいのではないかなと思いました。

LT内容

Google Presentationと事前に録音したYoutubeをこちらに載せておきます。

発表音声の文字起こし

LTを7分以内に抑えるリハーサルをした際の音声があったので、それをGoogle Speech To Text APIを使って文字起こしをしておきます。
「それではその音に近づくための研究と練習ということで江口さん自己紹介ですが名前は救えバックグランドは経済学や統計学仕事では同席部のマネージャー兼データ分析さんはこんな感じで趣味はトランペットブログ作成をやってるので良かったみてください結局どうしたんですけどその他データ分析関連のかって言いよったでどうやって勝手の音消したいなって思いますがここまでのホスピスがまずあって写真撮って進度が速い遅いである時と言うとって治ってるということですが音の長さが変わりますなのでお体変わります秘密のピストンだけなので発送の仕方ありませんのでコントロールして40種類ぐらいのレベルからは死んだことに帰って前記処理学会のそれの大会でスペクトルグラムを使ってプロのアクセスをしたのとの比較ではその方向の出来不出来を判定すると犬が行われていました鬼時間に周波数色の濃淡で信号ってなっちゃったから明日グラフとなりますがこの場合は時間が12秒まであって周波数が8000まで行ってあのね精神旺盛な強さが違ってます本当の鳴き声何もとしてるの俺も行ってここだったから体操マット耐熱耐酸プルマンちょっと寄りますよねプレイヤーの力も貸したして自分との違いを明らかにして違いがある場合とならそこによるそれが産めるのかお話していきたいと思ってます結果こんな感じえーとまずはサイトメトロ ffmpeg をインストールしてああいパイソン mrtg をインストールしたギターの音は愛想の良いしらすを呼び出してこういう形で書くことでノートブックので他のところでエロ金ボタンが現れますこれ音声は録音開始でもう一度録音するようになっていますその後に録音したファイルを保存 SMS でまあ16000hz 2のファイルに変換してまた再読み込みという形があってもはい録音環境としてはわからないぞマイクでして何楽器はちょっと古いので今回あったのがという音楽などの分析見ております玄関のために必要なセッティングってのがあるんでは基本的なセッティングで言っても実際に自分のそれかしたら私に出てきますね自分のやってみようと思いますまず聞いてみましょう第4フォルマント信号の強さがえっと Pro は赤いんですけど強いんですけど出してもらうわいえ私なってます遅く音源第4フォルマントの方がやっぱり行けないなっていうのが分かれば鼻が詰まってるから行けないっていましたねと自分ではまあ剛てる所間違いが分かりましたどうしたら強い信号になってるのかその女お腹に力を入れる行きの数量変えるとかありますもうなんかの力を入れて強いわ強い力を込めてもっと明るくならなかったにやと言われるとも出てこなくなったよ沢山入ったとこですねこれで第4子が赤くなるのはこれまで行きの量が足りてなかった自分じゃなかったじゃないかっていうので行ってみましょう早くなりましたねでスタグラムに従ってくると自分の音が帰宅したところ息を吸う量を多くすることで気分良くなることがわかりました待機するのって管楽器重要な今なってどれぐらい剃ったらいいのか分かんなかったり根性論だったりするんでこれは便利だねまあは管楽器にも分からなくなって鳴るんですけど G 計り方とかあると結婚もしないからっても行きたいのでくれてありがとうございました」

ソースコード（LT）

こちらではLTのために使ったスペクトログラムを表示するためのコードを共有します。

ソースコード（文字起こし）

Google Speech To Text APIで文字起こしをする際のコードをここに載せておきます。

GCSのバケットを作っておきます。今回はperry_voiceという名前のバケットを作っておきました。先ほど作ったaudio_only_filtered_441.wavを、このバケットにあげておきます。

APIの認証用のjsonファイルをcloud shellにアップロードしておきます。jsonファイルはGCPのコンソールの認証情報からサービスアカウントキーを作成し、「鍵を追加」の項目からjsonファイルを選択して手に入れることができます。（当然、Google Speech To Text APIのAPIを有効にしておく必要があります。）
1ヶ月につき60分までは無料なので、今回の処理も無料でできます。（GCSは数円かかると思いますが。）

文字起こし結果はGCSのperry_voiceの中にあります。

他の発表から学んだこと

深層学習を用いたJazz音楽の研究や、音楽の理論、DJプレイの自動化の研究、声の変換技術など同じ音という題材ながらも幅広く話を聞くことができました。私だけ入門な内容で申し訳ない気持ちが生じましたが、これからこの分野を深めるのもプレイヤーとして面白いなと思ったので、長く学ばせてもらおうと思います。

懇親会で学んだこと

・NVIDIAのノイズリダクション技術は文字起こしの際にやってみる価値あり
・賃貸の場合、防音の部屋は家賃相場の2~3万円アップ
・口周辺の筋肉の研究をされている方がいる（論文でその分野があるのは知っていたが、実験をやっていた方に会えた）
・みな、条件を揃えるためにマイクや集音の仕方にこだわりがあった。
・いい音の定義にやはり答えはない。

本来、正月の企画物としてアイデアを出していたものの、作成する時間がとれなかったのですが、日常的にもっと確率に従って生活できたら面白いなと思って作ってみることにしました。

二項分布

• やりたいこと
  オフィスにおいて、チームごとにデスク周辺の掃除当番を毎日乱数で決めているとして、1ヶ月で自分がどれくらい掃除をやらなくてはいけないのかを知りたい。あるいは、1回の掃除で済む確率とか、10回掃除する確率とかを知りたい。

  これは、成功確率p（掃除当番になる確率 = 1÷チームの人数 ）でのN回の独立なベルヌーイ試行となり、自分が掃除当番になった場合X=1で、そうでない場合はX=0となる。
  N回ベルヌーイ試行した際のXの合計は二項分布従うので、N回試行して、Xがkになる確率は以下のように示される。（チームの人数はHとする。）

  $$ P(Y_N = k ) = {}_N C _k \left (\frac{1}{H} \right )^k \left ( 1- \frac{1}{H} \right )^{N-k} \\
  (0 \leq k \leq N )$$
• コード
  
  from scipy.stats import binom import matplotlib.pyplot as plt H = 10 # チーム人数 D = 20 # 営業日 n, p = D, 1/H # 平均、分散、歪度、尖度 mean, var, skew, kurt = binom.stats(n, p, moments='mvsk') print('平均：',mean) fig, ax = plt.subplots(1, 1) x = range(D + 1) ax.plot(x, binom.pmf(x, n, p, loc=0), 'bo', ms=8, label='binom pmf') ax.vlines(x, 0, binom.pmf(x, n, p), colors='b', lw=5, alpha=0.5); ax.set_xlabel("count"); ax.set_ylabel("probability"); ax.grid(color='b', linestyle='-', linewidth=0.1); ax.set_title("Probability mass function");

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

  from scipy.stats import binom import matplotlib.pyplot as plt   H = 10 # チーム人数 D = 20 # 営業日   n, p = D,  1/H   # 平均、分散、歪度、尖度 mean, var, skew, kurt = binom.stats(n, p, moments='mvsk')   print('平均：',mean)   fig, ax = plt.subplots(1, 1)   x = range(D + 1) ax.plot(x, binom.pmf(x, n, p, loc=0), 'bo', ms=8, label='binom pmf') ax.vlines(x, 0, binom.pmf(x, n, p), colors='b', lw=5, alpha=0.5); ax.set_xlabel("count"); ax.set_ylabel("probability"); ax.grid(color='b', linestyle='-', linewidth=0.1); ax.set_title("Probability mass function");

  
  
  平均では20営業日で2回は掃除することになりますが、4回以上掃除しないといけない確率が10%以上はあることがわかります。掃除が嫌いな人は、苦しみの上限をイメージできると楽なのではないかと思うので、確率に従い掃除当番を甘んじて受け入れましょう。
• 参考情報
  scipy.stats.binom

matplotlib.pyplot.grid

多項分布

• やりたいこと
  先ほどのように乱数で掃除当番を割り当て、それをランキングにし、一位を掃除大臣として称号を付与するものとする。10名いるもののうち、上位のメンバーA・メンバーBとで累積の掃除回数がそれぞれ10回、6回としたもとで次の月に掃除大臣が入れ替わる確率を計算したい。
• コード
  
  import math H = 10 # チーム人数 D = 20 # 営業日 # あるメンバーがa回当たり、あるメンバーがb回当たる確率 def calc_prob_multinomial(a, b, H=H, D=D): return math.factorial(D) / (math.factorial(a) * math.factorial(b) * math.factorial(D - a - b)) * (1 / H)**(a) * (1 / H)**(b) * (1 - 2 / H )**(D - a - b) base_a = 10 # メンバーaの現在の掃除回数 base_b = 6 # メンバーbの現在の掃除回数 gap = base_a - base_b # 差分 prob_sum = 0.0 for i in range(D + 1): for j in range(1, D + 1): b = gap + i + j # 差分に1を足すことで負けるケースを計算（b = jにしたら確率の合計は1になった。） if b + i < D + 1: # 合計が月数を超えないようにする。 prob_sum = prob_sum + calc_prob_multinomial(i, b, H=H, D=D) # 確率を足す。 print('sum:',prob_sum)

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

  import math   H = 10 # チーム人数 D = 20 # 営業日   # あるメンバーがa回当たり、あるメンバーがb回当たる確率 def calc_prob_multinomial(a, b, H=H, D=D):     return math.factorial(D) / (math.factorial(a) * math.factorial(b) * math.factorial(D - a - b)) * (1 / H)**(a) * (1 / H)**(b) * (1 - 2 / H )**(D - a - b)   base_a = 10 # メンバーaの現在の掃除回数 base_b = 6 # メンバーbの現在の掃除回数 gap = base_a - base_b # 差分   prob_sum = 0.0 for i in range(D + 1):     for j in range(1, D + 1):         b = gap + i + j  # 差分に1を足すことで負けるケースを計算（b = jにしたら確率の合計は1になった。）         if b + i  < D + 1: # 合計が月数を超えないようにする。             prob_sum = prob_sum + calc_prob_multinomial(i, b, H=H, D=D) # 確率を足す。          print('sum:',prob_sum)

  
  計算したところ、1.2%くらいになった。このケースにおいてAさんがBさんに次の月に掃除大臣の座を奪われる可能性は低いらしい。
• 参考情報
  scipy.special.comb

誕生日の問題

• やりたいこと
  あまりにも有名な確率の話で、誕生日のパラドックスとも呼ばれています。誕生日が同じ人が部署内にいる確率が0.5を超える人数は何人か？
• コード
  
  import matplotlib.pyplot as plt N = 50 # メンバーの数の上限 D = 365 # 誕生日の数 def calc_prob_birth(N=N, D=D): prob = 1.0 for i in range(1, N): prob = (1 - i /D) * prob return 1 - prob fig, ax = plt.subplots(1, 1) x = range(2 ,N + 1) probability = list(map(lambda x:calc_prob_birth(N=x) , x)) ax.plot(x, probability, 'bo', ms=8) ax.hlines(y=0.5, xmin= 0, xmax=N, color='r'); ax.set_xlabel("the number of classmate"); ax.set_ylabel("probability"); ax.grid(color='b', linestyle='-', linewidth=0.1); ax.set_title("Probability");

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

  import matplotlib.pyplot as plt   N = 50 # メンバーの数の上限 D = 365 # 誕生日の数   def calc_prob_birth(N=N, D=D):     prob = 1.0     for i in range(1, N):         prob = (1 - i /D) * prob       return 1 - prob   fig, ax = plt.subplots(1, 1)   x = range(2 ,N + 1) probability = list(map(lambda x:calc_prob_birth(N=x) , x)) ax.plot(x, probability, 'bo', ms=8) ax.hlines(y=0.5, xmin= 0, xmax=N, color='r'); ax.set_xlabel("the number of classmate"); ax.set_ylabel("probability"); ax.grid(color='b', linestyle='-', linewidth=0.1); ax.set_title("Probability");

  
  
• 参考情報
  誕生日のパラドックス
  matplotlib.pyplot.hlines

ファーストサクセス分布

• やりたいこと
  アルバイトの面接を受けるとして、その合格率が1/30だとして、何回くらい面接しても合格できないものだろうか。苦しみの上限を知りたい。20回、30回面接を受けた時に合格できてない確率はどれくらいだろうか。何回くらい面接を受けたら合格できていない確率を1%未満に下げれるだろうか。
  成功確率がpのベルヌーイ試行を独立に何回も行うとき、初めて成功数までに失敗した回数をXとして、初めて成功した際の試行回数をYとした場合、その確率は\( k \geq 1 \)として、
  $$ P(Y = k) = P(X=k-1)=p(1-p)^{k-1} $$

  と表すことができる。これをパラメータpのファーストサクセス分布と呼ぶ。（幾何分布とも呼ぶ。）

• コード
  J = 1/ 30 # 面接合格率 def calc_prob_first_success(J=J, k=k): return (J) * (1- J )**(k-1) fig, ax = plt.subplots(1, 1) x = range(1 , 10000) # 無視できるくらい伸ばしておく。 probability = list(map(lambda x:calc_prob_first_success(k=x,) , x)) probability_sum = list(map(lambda x:1 - sum(probability[:x]) , x)) length = 90 # 90回の失敗までを見る。 ax.plot(x[:length], probability_sum[:length], 'bo', ms=8) ax.set_xlabel("the number of failure"); ax.set_ylabel("probability"); ax.grid(color='b', linestyle='-', linewidth=0.1); ax.set_title("the probability we can not success on the number of failure");

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

  J = 1/ 30 # 面接合格率   def calc_prob_first_success(J=J, k=k):     return (J) * (1- J )**(k-1)   fig, ax = plt.subplots(1, 1)   x = range(1 , 10000) # 無視できるくらい伸ばしておく。 probability = list(map(lambda x:calc_prob_first_success(k=x,) , x)) probability_sum = list(map(lambda x:1 - sum(probability[:x]) , x)) length = 90 # 90回の失敗までを見る。   ax.plot(x[:length], probability_sum[:length], 'bo', ms=8) ax.set_xlabel("the number of failure"); ax.set_ylabel("probability"); ax.grid(color='b', linestyle='-', linewidth=0.1); ax.set_title("the probability we can not success on the number of failure");

  

  20回面接をしても合格できない確率はおよそ50%。30回面接しても合格できない確率はおよそ36%。合格できない確率が1%を切るのは面接回数が135回になった時。このレベルを追い求めると心が折れる面接回数なのかもしれない。

• 参考情報
  弱点克服 大学生の確率・統計

負の二項分布

• やりたいこと
  ある営業がリンゴを売る仕事をしておりタウンページで飲食店に電話をかけているとする、商談の予約が取れる確率(p)が0.05として、商談予約が5(n)回起きるまでに失敗がk回生じる確率について知りたい。いったい、この営業は飲食店にどれくらい商談をお断りをされる覚悟で臨めばいいのか。

  このような事象は負の二項分布に従うので、
  $$ P(X = k) = {}_{n – 1 + k} C _k p^{n-1}(1-p)^k $$
  という確率に従う。
  なお、この分布の期待値は
  $$ E(X) = \frac{n(1-p)}{p} $$
  に従う。

• コード
  import matplotlib.pyplot as plt import math prob = 1/ 20 # リンゴ営業の商談予約できる確率 n = 5 # 営業の目標である商談予約の回数 # 負の二項分布の確率密度関数 def calc_prob_nb(prob=prob, K=K, n=n): return math.factorial(K+n-1) / ( math.factorial(n-1) * math.factorial( K+n-1 - (n-1) )) *(prob)**n * ( 1- prob )**K fig, ax = plt.subplots(1, 1) x = range(1 , 200) probability = list(map(lambda x:calc_prob_nb(prob=prob, K=x, n=n) , x)) length =200 ax.plot(x[:length], probability[:length], 'bo', ms=8) ax.set_xlabel("the number of failure"); ax.set_ylabel("probability"); ax.grid(color='b', linestyle='-', linewidth=0.1); ax.set_title("Probability until the success happens on five times");

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

  import matplotlib.pyplot as plt import math   prob = 1/ 20 # リンゴ営業の商談予約できる確率 n = 5 # 営業の目標である商談予約の回数   # 負の二項分布の確率密度関数 def calc_prob_nb(prob=prob, K=K, n=n):     return math.factorial(K+n-1) / (  math.factorial(n-1) * math.factorial(  K+n-1 - (n-1) )) *(prob)**n * ( 1- prob )**K   fig, ax = plt.subplots(1, 1)   x = range(1 , 200) probability = list(map(lambda x:calc_prob_nb(prob=prob, K=x, n=n) , x)) length =200   ax.plot(x[:length], probability[:length], 'bo', ms=8) ax.set_xlabel("the number of failure"); ax.set_ylabel("probability"); ax.grid(color='b', linestyle='-', linewidth=0.1); ax.set_title("Probability until the success happens on five times");

  
  75回くらいの失敗で確率自体のピークがきている。期待値だと95回だが、確率のピークはもっと手前にある。

  累積確率を知りたい場合はこちらのコード

  fig, ax = plt.subplots(1, 1) x = range(1 , 200) # 無視できるくらい伸ばしておく。 probability = list(map(lambda x:calc_prob_nb(prob=prob, K=x, n=n) , x)) probability_sum = list(map(lambda x:sum(probability[:x]) , x)) length =200 ax.plot(x[:length], probability_sum[:length], 'bo', ms=8) ax.set_xlabel("the number of failure"); ax.set_ylabel("probability"); ax.grid(color='b', linestyle='-', linewidth=0.1); ax.set_title("Cumulative probability until the success happens on five times");

  1 2 3 4 5 6 7 8 9 10 11

  fig, ax = plt.subplots(1, 1)   x = range(1 , 200) # 無視できるくらい伸ばしておく。 probability = list(map(lambda x:calc_prob_nb(prob=prob, K=x, n=n) , x)) probability_sum = list(map(lambda x:sum(probability[:x]) , x)) length =200   ax.plot(x[:length], probability_sum[:length], 'bo', ms=8) ax.set_xlabel("the number of failure"); ax.set_ylabel("probability"); ax.grid(color='b', linestyle='-', linewidth=0.1); ax.set_title("Cumulative probability until the success happens on five times");

  
  80%のところで125回なので、5人に一人は期待値の95回を30回くらいは超えてしまってもおかしくないということになる。そして、5人に一人くらいは50〜60回で達成できるラッキーな人もいる。

• 参考情報
  弱点克服 大学生の確率・統計

ポアソン分布

• やりたいこと
  先ほど、営業が飲食店に売ったリンゴの中に芋虫が住み付いているケースが500回に1回の頻度で発生し、営業と顧客の間でトラブルに発展するとする。年間に1200の飲食店と取引を行うが、トラブルはどれくらい発生するものなのか。
  顧客とのトラブルが発生する回数がYだとした際に、ポアソン分布に従うとした場合、
  $$P(Y = k) = \frac{\lambda^k}{k!} e^{-\lambda} $$
  となる。ここで、
  $$ \lambda = 1200 \times \frac{1}{500} $$
  とおく。
  なお、期待値と分散はともに
  $$ E(Y) = V(Y) = \lambda $$
  となる。
• コード
  
  import matplotlib.pyplot as plt import math import numpy as np lamda = 1200 * (1 /500 ) def calc_prob_poisson(lamda=lamda, K=K): return lamda**K / math.factorial(K) * np.exp(-lamda) fig, ax = plt.subplots(1, 1) x = range(0 , 10) probability = list(map(lambda x:calc_prob_poisson(lamda=lamda, K=x) , x)) length =10 ax.plot(x[:length], probability[:length], 'bo', ms=8) ax.set_xlabel("the number of trouble"); ax.set_ylabel("probability"); ax.grid(color='b', linestyle='-', linewidth=0.1);

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

  import matplotlib.pyplot as plt import math import numpy as np   lamda = 1200 * (1 /500 )   def calc_prob_poisson(lamda=lamda, K=K):     return lamda**K / math.factorial(K) * np.exp(-lamda)   fig, ax = plt.subplots(1, 1)   x = range(0 , 10) probability = list(map(lambda x:calc_prob_poisson(lamda=lamda, K=x) , x)) length =10   ax.plot(x[:length], probability[:length], 'bo', ms=8) ax.set_xlabel("the number of trouble"); ax.set_ylabel("probability"); ax.grid(color='b', linestyle='-', linewidth=0.1);

  
  
  トラブルが2回の時が確率が一番高い。0回で済む確率も0.1に近いくらいはある。
  
  fig, ax = plt.subplots(1, 1) x = range(0 , 10) probability = list(map(lambda x:calc_prob_poisson(lamda=lamda, K=x) , x)) probability_sum = list(map(lambda x:sum(probability[:x]) , x)) length = 10 ax.plot(x[:length], probability_sum[:length], 'bo', ms=8) ax.set_xlabel("the number of trouble"); ax.set_ylabel("probability"); ax.grid(color='b', linestyle='-', linewidth=0.1); ax.set_title("Cumulative probability");

  1 2 3 4 5 6 7 8 9 10 11

  fig, ax = plt.subplots(1, 1)   x = range(0 , 10) probability = list(map(lambda x:calc_prob_poisson(lamda=lamda, K=x) , x)) probability_sum = list(map(lambda x:sum(probability[:x]) , x)) length = 10   ax.plot(x[:length], probability_sum[:length], 'bo', ms=8) ax.set_xlabel("the number of trouble"); ax.set_ylabel("probability"); ax.grid(color='b', linestyle='-', linewidth=0.1); ax.set_title("Cumulative probability");

  
  
  5人に一人は4回以上のトラブルに巻き込まれる気の毒な営業がいそうですね。まぁ、リンゴの中は開けてみないとわかりませんから、顧客トラブルがあっても恨みっこなしで頑張りましょう。
• 参考情報
  弱点克服 大学生の確率・統計

フィギュア集め問題

• やりたいこと
  現場猫フィギュアのガチャポンに関心があるとする。N回ガチャポンを回した時に、N種類の現場猫フィギュアが全部揃う確率はどれくらいだろうか？（どの現場猫のフィギュアも出現確率は同じものとする。）
  また、N+1回ガチャポンを回したとき、N種類の全部が揃う確率はどれくらいだろうか？あるいは、m回ガチャポンを回してそろったフィギュアの数の期待値は？そして、Y回目のガチャポンではじめてN種類の現場猫の全部揃う、Yの期待値は？

  ・N回ガチャポンを回した時に、N種類の現場猫フィギュアが全部揃う確率：
  $$\frac{N!}{N^N}$$

  ・N+1回ガチャポンを回したとき、N種類の全部が揃う確率：
  $$ N \times \frac{(N+1)!}{2} \times \frac{1}{N^{N+1}} $$

  ・m回ガチャポンを回してそろったフィギュアの種類\( X_m \)の期待値：
  m回ガチャポンを回してそろったフィギュアの種類\( X_m \)を以下のように表す。
  $$X_m = Z_1 + Z_2 + \dots + Z_N$$
  m回ガチャポンを引いて、当たらないということは
  $$ P(Z_i = 0) = \left ( \frac{N – 1}{N} \right )^m $$
  ということになる。
  これをもとに\( Z_i \)の期待値を計算すると

  $$ E(Z_i) = 1 \times \left ( 1 – \left ( \frac{N – 1}{N} \right )^m \right ) + 0 \times \left ( \frac{N – 1}{N} \right )^m \\
  = 1 – \left( \frac{N-1}{N} \right )^m $$
  よって、\( X_m \)の期待値は
  $$ E(X_m) = E( Z_1 + \dots + Z_N) \\
  = N \left ( 1 – \left( \frac{N-1}{N} \right )^m \right )$$

  となる。

  ・Y回目のガチャポンではじめてN種類のフィギュアが揃う際のYの期待値：
  i-1種類からi種類になるまでのガチャポンの回す回数を\( T_i \)として、その階差がファーストサクセス分布に従うとする。
  $$ T_i – T_{i-1} \sim F_s \left ( \frac{N – i + 1}{N} \right ) $$
  その場合、Yの期待値は以下のようになる。

  $$ E(Y) = E( T_N ) \\
  = E( (T_{N} – T_{N-1}) + (T_{N-1} – T_{N-2}) + \dots + (T_{2} – T_{1}) + T_{1} ) \\
  = E \left ( F_s \left ( \frac{1}{N} \right ) + F_s \left ( \frac{2}{N} \right ) \dots + F_s \left ( \frac{N-1}{N} \right ) + 1 \right ) \\
  = N \left ( 1 + \frac{1}{2} + \dots + \frac{1}{N} \right )
  $$
• コード
  
  import matplotlib.pyplot as plt import math N = 5 # 現場猫フィギュアの種類 m = 5 # 回すガチャポンの回数m # N回ガチャポンを回した時に、N種類の現場猫フィギュアが全部揃う確率 print(math.factorial(N) / (N ** N)) # N+1回ガチャポンを回したとき、N種類の全部が揃う確率 print(N * math.factorial(N+1) / 2 * (1/N**(N+1)) ) # m回ガチャポンを回してそろったフィギュアの数のXの期待値 print(N * (1 - ((N-1)/N)**m )) # Y回目のガチャポンではじめてN種類のフィギュアが揃う際のYの期待値 def harmonic_series_sum(n): i = 1 s = 0.0 for i in range(1, n+1): s = s + 1/i; return s; print(N * (harmonic_series_sum(N)))

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

  import matplotlib.pyplot as plt import math   N = 5 # 現場猫フィギュアの種類 m = 5 # 回すガチャポンの回数m   # N回ガチャポンを回した時に、N種類の現場猫フィギュアが全部揃う確率 print(math.factorial(N) / (N ** N))     # N+1回ガチャポンを回したとき、N種類の全部が揃う確率 print(N * math.factorial(N+1) / 2 * (1/N**(N+1)) )     # m回ガチャポンを回してそろったフィギュアの数のXの期待値 print(N * (1 - ((N-1)/N)**m  ))     # Y回目のガチャポンではじめてN種類のフィギュアが揃う際のYの期待値 def harmonic_series_sum(n):     i = 1     s = 0.0     for i in range(1, n+1):         s = s + 1/i;     return s;   print(N * (harmonic_series_sum(N)))

  
  N回ガチャポンを回した時に、N種類の現場猫フィギュアが全部揃う確率：
  0.0384
  N+1回ガチャポンを回した時、N種類の現場猫フィギュア全部が揃う確率：
  0.1152
  m回ガチャポンを回して揃った現場猫フィギュアの種類Xの期待値：
  3.3615
  Y回目のガチャポンではじめてN種類の現場猫フィギュアが揃う際のYの期待値：
  11.4166
• 参考情報
  弱点克服 大学生の確率・統計

以前、トランペットのマウスピース選びで音声解析をしてみましたが、今回は楽器選びで行ってみます。
前回の調査でf1値とかでは良い音の判断がつかなかったので、スペクトログラムを出してみたいと思います。
管楽器はそこそこに高い買い物なので、色々と比較見当をしてから意思決定したいですね。

今回の楽器たちはヤマハの800GS（所持品）、Bach 180ML37SP、Bach 180ML37GL、Bach AL190GLです。ヤマハのやつは20万円くらいですが、Bachのは30万円くらいします。ちなみにBachはアメリカの楽器メーカーで、私の感覚値ではあるのですが、音大生はだいたい持っているような気がします。10万円の差は大きいのか小さいのか。そしてプロの演奏家のスペクトグラムとも比較します。

スペクトログラムとは

パワースペクトルの値の高低を画像の濃淡で表現し、縦軸に周波数、横軸に時間を取った座標上にプロットしたものをスペクトグラムと呼びます。スペクトグラムでは、周波数が低いところ（下）から、第1フォルマント、第2フォルマント、第3フォルマント・・・と呼びます。フォルマントは周りの周波数と比べて突出している部分のことを指します。

$$ S(\omega) = \int_{-\infty}^{\infty}s(t) e^{-jft}dt$$

さて、いきなりの数式ですが、フーリエ変換の式となります。\( S(\omega) \)を二乗したものがパワースペクトルとなります。（二乗するは複素平面上で絶対値を取りたいため。）
tは時刻、fは基本周波数（最も低い周波数）、jは虚数単位を表します。ここで登場する\( s(t) \)は定常波と呼ばれ、以下の式で定まります。

$$s(t) = \frac{a_0}{2} + \sum_{n=1}^{\infty}a_n \cos (2 \pi n ft + \theta_n)$$
\( a_n \)は元の信号\( s(t) \)にどの周波数の波がどの程度の強さで含まれているかという情報を表します。
フーリエ変換のモチベーションとしては、周波数の周期関数を周波数の異なる単純な波の重み付き和で表現できることですが、トランペットの音波を色々な周波数に分解できるのは面白いです。良い奏者の周波数がどうなっているのか調べてみるのは面白そう。

Pythonでの計算&比較

音声データはFFMpegで変換とか切り取りをするなどして作成しました。作成方法は簡単で、以前のブログにも記しましたので、試してみたい方はそちらを参考にしてみてください。

まずはヤマハですが、この楽器になれている分、スムーズに吹けています。最初のほうのリップスラーは5つくらいのフォルマントがありそうです。音階が上がるにつれてフォルマントが減っているのがわかります。最後のHigh Bでは細めのフォルマントが3つあります。


次は、Bach 180ML37SPという楽器です。鳴らすのに慣れていないので、ゆっくり吹きました。最初のほうのリップスラーは5つくらいのフォルマントがありそうで、それはヤマハと同じですが、最後のHigh Bで2つくらいしかフォルマントがありません。楽器に十分息が吹き込めておらず、鳴らせていないのだろうと思われます。可視化できるって面白いですね。


次は、Bach 180ML37GLという楽器です。さっきの楽器との違いは銀メッキがなされているかどうかくらいです。最初のリップスラーは同じ感じで5つくらいのフォルマントがありそうですが、リップスラーで戻ったところで5つ目のフォルマントが消えています。ただ、最後のHigh Bで3つほどフォルマントがありそうです。慣れてきたのでしょうか？


次は、Bach AL190GLという高級モデルです。40万円くらいします。より技術力の高い職人が手掛けた楽器だそうです。ベルに一手間加わっていたり、専用パーツがあったりするみたいです。実際、30万円のよりも吹きやすい印象がありました。
最初のリップスラーは同じ感じで5つくらいのフォルマントがありそうですが、最後のHigh Bでフォルマントが2つになっています。やっぱ鳴らすの難しいのでしょうか。


最後に、プロの方の音源です。最初のBで私のグラフと明らかにフォルマントの数が違うことがわかります。楽器を鳴らせているというのはこういうレベルなのかなと思ってしまいました。High Bでフォルマントが4つはありそうです。このレベルになるにはどういう吹き方をすればいいのか、色々研究のしがいがありそう。

参考文献

イラストで学ぶ 音声認識 (KS情報科学専門書)