遊び tokidoki 仕事

数学と音楽と教育と遊び

Top | おしごと | ゼミ | がくせい | すうがく | かがく | きょういく | おんがく | おきにー | Tips | Photo | イベント | ものもう | あれこれ | About

Benford則を見てみた

今年も講義「統計とコンピュータ」がはじまった.
それが一般の学生なら興味を持つ話題であっても数学の学生はデータを眺めることが嫌いだ,というこの十数年の経験が気持ちを曇らせるのだけど,また新たなネタ作りを試みている.
今年はBenford則をやってみよう.というより統計教育でよく取り扱われる標準ネタであるものを今までやっていなかったことに気付いた.
ベンフォードの法則 - Wikipedia

下準備で比較的大きな桁数まで現れそうな株式市場からのデータに適用してみた.
2018年3月30日の東証一部上場企業4019社の「始値」と「出来高」それぞれの最高桁の数字の分布を眺める.
始値については「株価」としての適正な価格範囲が共通認識としてあるために理論的になるはずの分布からややずれる一方で,出来高については発行株式数が巨大な企業もあるために,Benford則が綺麗に表れている.
おお,ネタにするには良い感じだ.
f:id:okiraku894:20180413232335p:plain

一方,数学の学生が興味を持ちそうな,与えられた数のべき乗の先頭の数字についてもやってみた.
(1.1)^n から (1.5)^nn=1 から 1000 まで行って各々先頭の数字を拾った.
こちらは更に見事にBenford則に当てはまっている.
f:id:okiraku894:20180413232339p:plain

結局のところ,これらの現象は確率変数  X\log_{10}X の小数部分が [0,1) で一様分布することが元になっているのであって,何か珍しいことが起こっているわけではないのだけど,見た目にはへ~っとなりやすい.

そもそも数  X の最高桁の数が  m だということは適当な自然数 e によって
\begin{equation}
m\cdot 10^e\le X < (m+1)\cdot 10^e
\end{equation}と書けるということであり,
\begin{equation}
e+\log_{10}m \le \log_{10}X < e+\log_{10}(m+1)
\end{equation}となるから,つまり \log_{10}X の小数部分を眺めているに他ならない.そして
\begin{equation}
\sum_{m=1}^9\left(\log_{10}(m+1)-\log_{10}m\right)=1
\end{equation}となるものだから,ちょうど  \log_{10}(m+1)-\log_{10}m が数字  m の現れる確率となる.

特に数  \alpha\ (\log_{10}\alpha\not\in \mathbb{Q}) のべき乗ならば, n\log_{10}\alpha [0,1) に一様分布するという良く知られたWeylの定理があることから,殊更Benford則がはっきり表れることになる.

さて,このネタ,反応あるかなぁ...

ようやく卒論審査会 in 2018

今年度もようやくここまで来ました.
f:id:okiraku894:20180112130147j:plain:w400,left


f:id:okiraku894:20180112155016j:plain:w400,right
タイトル,一人だけまだ決まってないとき.


緊張の本番開始.
f:id:okiraku894:20180210092826j:plain:w600

f:id:okiraku894:20180210093855j:plain:w400,left
結構うすい環(わ)ぁ,土星の環
―重力ポテンシャルを利用した惑星環形成のモデル―

土星などの惑星環が薄く形成されてしまう理由を簡単な力学だけで説明を試みたよ.


f:id:okiraku894:20180210094542j:plain:w400,right
タネト――ク
―カードマジックを離散力学系で読み解く―

カードシャッフル・トリックにまつわる数理的現象を離散力学系として分析した.その場でマジックをする初めての試み.


f:id:okiraku894:20180210100423j:plain:w400,left
matrix( 母 ) をたずねて三次元
―高次元ピタゴラス数のネットワーク―

ピタゴラス数を産む有名な行列の話を3次元にも拡張できないか試みた話.今回は3次元にお母さんは見つからなかった.果たして3次元以上にお母さんはいるのだろうか?


f:id:okiraku894:20180210101341j:plain:w400,right
効率ファースト
―待ち行列を用いたより効率的なレジ待ち並びの提案―

待ち行列理論をレジ待ちに適用した.あんなに苦労した不等式評価がさらっと登場しただけなんだけど,発表するとやはりそうなるのかな.


f:id:okiraku894:20180210103112j:plain:w400,left
ゼミ友学園
―研究室マッチング問題,忖度を添えて―

色んな意味でブラックジョークを交えて,1対多マッチング理論とその応用.Scratchによるアルゴリズムの提示を今回初めてやってみた.


f:id:okiraku894:20180210105040j:plain:w400,right
ゼミの中心で「あ~」と叫ぶ
―線形予測を用いた sin 波の合成―

Excel VBAも使わず,ほとんど手作業で線形予測による母音の合成を行うという力作.いやホント,合成が間に合ってよかったよ.ちゃんと「あいうえお」って聞こえたし,ね.


うん,何とかなった,今年も.皆の者,乙.



他のゼミの発表を見て改めて考えさせられる点,多々あり.
徐々にうちのゼミも下学年に「分かった気にさせる」方向を取ってきたのだが,
果たしてそれをどこまですべきなのか,ということだった.
つまりこれら発表が「誰」あるいは「何」に対して行われるべきか,ということだ.
他ならぬこれは「卒業論文審査会」であって,ゼミの宣伝会ではない.
審査員向けの研究発表が本題であることに変わりはなく,
その点に真摯に向き合っているかどうか,は守らねばならない点だ.

しかし一方で下の学年の勉学への意欲向上につなげる,という目的もこの審査会にはある.
(と思ってたけど,そう?)
分かりやすいところ,受けの良いところだけを掻い摘んで発表するスタイルも,
真の研究内容の誤解を招かない範囲で行われても可,とも思える.
分からない話が延々と続くのであれば,彼らの意識も持たない.
ちょいちょい下学年向けの小ネタを挟みながらのトークも必要悪なのかもしれない.
さて,ではどうあるべきだろうか?

下学年には悪いのだけど,表面だけ「分かったつもり」になってもらうこと目指しながらも,
(それはきっと十分な間をかけて研究動機を語ることなのだろう)
分かっている人が見れば研究内容の本質がちゃんと語られている,
そんなプレゼンが落としどころなんだろうかね,という気になってきた.
この辺りのバランスが難しいし,教員目線だけでは分からないところだったりする.
ああ,でもそうだな,自分だったらこのネタ,どうプレゼンするかな,って思えば行ける気がする.

どちらにしてもそろそろBeamerを超えた,視覚に訴えるプレゼン道具が必要に思えた.
(え,なに,結局keynoteが羨ましい,という話かい?)


って,何とか終わって,う・ち・あ・げ!
でも,2年の世話人だったので開始から1時間遅れで参加.
そうそう,3年とそして去年卒業の11代目もゲスト参加.
f:id:okiraku894:20180212100226j:plain
f:id:okiraku894:20180212100225j:plain

で,口上と共に飲み会の場でいただいたのが,
若かりし頃(1年半前)の彼らの写真の入ったマグカップ,
f:id:okiraku894:20180212094846j:plain

から~の,ど~ん.

f:id:okiraku894:20180212095139j:plain:w400
ま,散々,九平次九平次ってゆうてたからね,ゼミ中に.
あざっす!

九平次とブルックナーと,そして卒論添削

今年の暮れも,例年の如くゼミ生の卒論添削で年を越す.
f:id:okiraku894:20171221114638j:plain
こうして手直しをする中,結構面白い発見をする.
例えばピタゴラス数を生成する行列の話は,
ピタゴラス数を2パラメータで表示するとき,そのパラメータが動く空間と
原始ピタゴラス数全体のなす空間の関係は,
パラメータ(2,1)を始対象とし,射を(1,1)(2,0)とする圏と
原始ピタゴラス数を含む圏との関手\Phi_2を観察してるに他ならないことだとか,
f:id:okiraku894:20171231202503p:plain
カードマジックで有名なGergonneトリックは一度論文にしたけれど,
Generalized Gergonne's Trick and its Continuous Approximation
いま一度卒論ネタにしてまとめ直すと綺麗に書き直せることとか.
f:id:okiraku894:20171231203127p:plain

まぁ,何にしても年末,
すき焼きとネットでわざわざ注文して手に入れた九平次(雄町!)の取り合わせを楽しんでから,
f:id:okiraku894:20171231175422j:plain:w500
昨年同様,遥か遠く高校時代に衝撃を受けた
マタチッチ指揮のN響によるブルックナー8番を聴きながら,卒論の添削で年を越すのだ.

ブルックナー:交響曲第8番

ブルックナー:交響曲第8番

  • アーティスト: ロヴロ・フォン・マタチッチ,NHK交響楽団
  • 出版社/メーカー: 日本コロムビア
  • 発売日: 2011/03/23
  • メディア: CD
  • 購入: 1人 クリック: 1回
  • この商品を含むブログを見る
悪くは無い.むしろ誰の為でもなくある種の芸術活動として没頭できるこの時間は
おそらく誰にも理解されないであろう,私だけの幸福のひと時である.

sin波による疑似フォルマントシミュレーション

今年の4年ゼミ生の一人は音声分析・合成をテーマにしている.
EXCEL上でちまちま作業してもらっているが,
線形予測によるフォルマント抽出がうまくいかないのか,
なかなか母音「あいうえお」すら \sin波から合成できないでいる.

手軽にBASICで母音の合成ぐらいできないものかと探したら,
高々6つ程度の倍音の重ね合わせで「あいうえお」らしい音を作っているページを見つけた.
qiita.com
本当に手軽にできそうだったので,10進BASICでサウンド生成する方法を探して作ってみた.
十進BASIC−外部プログラムの利用

でもって,音程も適当に定めて実行したのが↓

う~ん,そう思って聞けば聞こえなくもない程度.

AquesTalkつかえば,簡単なんだけどね↓

https://www.a-quest.com/index.html


以下,BASICソースと使い方.
実行すると,母音と音程の列を入力するよう聞いてくるので,
母音[aiueo]と音程[0~9]の組合せをinput.
上掲の音源では
   a3i7u9e0o2
とinputしてある.
するとソースのある場所に test.wav が生成されダイレクトにそれが再生される.
いや,その予定だったのだけど,やっているうちにダイレクト再生できなくなった.あれれ.
まぁ,生成された test.wav をクリックすれば再生できるようなので.

REM [疑似母音フォルマント シミュレータ]
REM Ver. 2017/12/23
REM 参考:
REM BASIC 元ネタ:http://hp.vector.co.jp/authors/VA008683/ExtProg.htm
REM               testwave.bas  十進BASICによる音声ファイル生成プログラム 
REM               chr$関数の引数として全ての1バイト整数が許される様に、オプションメニューの
REM               「互換性」「動作」で「文字列処理の単位」を「バイト」に設定して実行すること
REM フォルマント:https://qiita.com/rild/items/339c5c36f4c1ad8d4325

OPTION CHARACTER BYTE
DECLARE FUNCTION formant

INPUT  PROMPT "母音とその高さの列をどうぞ.(「あい↑う↑↑え↓お」は a3i7u9e0o2 など)":vw$
LET vwlen=LEN(vw$)
LET d=0.3  !継続時間(秒)
LET f0=220 !正弦波周波数(Hz)
LET db=-10 !生成波形の最大値(±32767を0dBとする)
LET fs=44100 !標本化周波数(Hz)
LET bps=fs*4 !1秒当りのデータ量(ステレオ16ビット量子化)
LET dsize=d*vwlen*fs*4 !オーディオデータサイズ
LET fsize=dsize+36 !ファイルサイズ(先頭8Bを除く)
LET fmtsize=16 !フォーマットサイズ
LET channel=2^17+1 !ステレオPCMデータの指定
LET reso=2^20+4 !16ビットの指定
LET a=10^(db/20)
OPEN #6 : NAME "test.wav" !ファイルを開き、この名前のファイルが
ERASE #6 !既に存在していた場合には上書きを指定
LET t0=TIME
PRINT #6 : "RIFF"; !以下、wavファイルのヘッダーを作成
CALL out4(fsize) !
PRINT #6 : "WAVEfmt "; !
CALL out4(fmtsize) !
CALL out4(channel) !
CALL out4(fs) !
CALL out4(bps) !
CALL out4(reso) !
PRINT #6 : "data"; !
CALL out4(dsize) !ここまでがヘッダー用の出力
LET al=a*32767 !左チャネル係数
LET ar=a*65535 !右チャネル係数
LET k0=f0/fs*PI*2 !引数の刻み
LET audio$="" !オーディオデータバッファを初期化
FOR vw=0 TO vwlen/2-1
   LET count=0 !カウンタをリセット
   LET vow$=vw$(vw*2+1:vw*2+1)
   LET ht=VAL(vw$(vw*2+2:vw*2+2))
   PRINT vow$,ht
   FOR i=1 TO d*fs
      LET vsin=formant(vow$,k0*(1+ht/20))
      LET lch=INT(vsin*al+0.5)
      IF lch<0 THEN LET lch=lch+65536 !負の数は補数表現
      LET rch=INT(vsin*ar+0.5)
      IF rch<0 THEN LET rch=rch+65536 !負の数は補数表現
      LET audio$=audio$&CHR$(MOD(lch,256))&CHR$(INT(lch/256))&CHR$(MOD(rch,256))&CHR$(INT(rch/256))
      LET count=count+1
      IF count=64 THEN !バッファデータが所定の長さ(64サンプル256Bがほぼ最適)に達したら
         PRINT #6 : audio$; !データをファイルに出力して
         LET audio$="" !バッファを初期化し
         LET count=0 !カウンタをリセット
      END IF 
   NEXT i
   PRINT #6 : audio$; !バッファに残ったデータを出力して
NEXT vw
CLOSE #6 !ファイルを閉じる
PRINT "elapsed time = ";TIME-t0;"seconds"

PLAYSOUND "test.wav"

SUB out4(i4) !4バイト整数の出力(little endian)
   LET j4=i4
   FOR m=0 TO 3
      PRINT #6 : CHR$(MOD(j4,256));
      LET j4=INT(j4/256)
   NEXT m
END SUB

FUNCTION formant(v$,k)
   SELECT CASE v$
   CASE "a"
      LET formant=0.19*SIN(i*k*4)+0.09*SIN(i*k*2)+0.08*SIN(i*k*3)+0.08*SIN(i*k*5)+0.07*SIN(i*k*1)+0.07*SIN(i*k*6)
   CASE "i"
      LET formant=0.19*SIN(i*k)+0.09*SIN(i*k*2)+0.08*SIN(i*k*11)+0.08*SIN(i*k*13)+0.07*SIN(i*k*12)
   CASE "u"
      LET formant=0.19*SIN(i*k)+0.09*SIN(i*k*6)+0.08*SIN(i*k*2)+0.08*SIN(i*k*5)+0.08*SIN(i*k*4)
   CASE "e"
      LET formant=0.19*SIN(i*k)+0.09*SIN(i*k*2)+0.08*SIN(i*k*3)+0.08*SIN(i*k*11)+0.07*SIN(i*k*4)
   CASE "o"
      LET formant=0.19*SIN(i*k*4)+0.09*SIN(i*k*2)+0.08*SIN(i*k)+0.08*SIN(i*k*3)
   CASE ELSE
      LET formant=SIN(i*k)
   END SELECT
END FUNCTION
END 

こんなところに Cayley-Hamilton

4年ゼミの話題の一つが,細矢氏が提案するトポロジカルインデックスなんだが,

トポロジカル・インデックス: フィボナッチ数からピタゴラスの三角形までをつなぐ新しい数学

トポロジカル・インデックス: フィボナッチ数からピタゴラスの三角形までをつなぐ新しい数学

その中にちょっと「おやっ」となった計算があったので,その紹介と種明かし.

ピタゴラス数を生成する行列にまつわる連立漸化式
\[
\begin{pmatrix}
a_{n+1}\\
b_{n+1}\\
c_{n+1}
\end{pmatrix}
=
A\begin{pmatrix}
a_n\\
b_n\\
c_n
\end{pmatrix}
=
\begin{pmatrix}
1 & 2 & 2\\
2 & 1 & 2\\
2 & 2 & 3
\end{pmatrix}
\begin{pmatrix}
a_n\\
b_n\\
c_n
\end{pmatrix}
\]から一種類の文字に関する漸化式
\[
a_{n+3}=5a_{n+2}+5a_{n+1}-a_n
\]等を導く件.本書では
\[
\sigma a_n=a_{n+1}
\]なる線形作用素を持ち出して
\[
\sigma\begin{pmatrix}
a_n\\
b_n\\
c_n
\end{pmatrix}
=A\begin{pmatrix}
a_n\\
b_n\\
c_n
\end{pmatrix},
\]すなわち
\[
(\sigma I-A)
\begin{pmatrix}
a_n\\
b_n\\
c_n
\end{pmatrix}
=\textbf{0}
\](Iは単位行列)がnontrivialな解をもつことから
\[
\det(\sigma I-A)=\sigma^3-5\sigma^2-5\sigma+1=0
\]が導かれるので
\[
(\sigma^3-5\sigma^2-5\sigma+1)a_n=a_{n+3}-5a_{n+2}-5a_{n+1}+a_n=0
\]なる漸化式が得られる,とある.
作用素\sigmaをあたかも数のように扱って議論しているわけだ.
はて,こんな扱いかた,妥当なんだろうか???



数日悩んだけど,何だか分かったので備忘録.
線形代数でいうところのCayley-Hamiltonの定理の証明を真似をすれば良い.
その証明は特性多項式を
\[
p(\lambda)=\det(\lambda I-A)
\]とし,(\lambda I-A)の余因子行列を\tilde{A}(\lambda)とおけば
\[
\tilde{A}(\lambda)(\lambda I-A)=p(\lambda)I
\]となるところから始まる.
さて,この両辺はともに\lambdaの多項式を成分とする行列なので,\lambdaへ作用素\sigmaを代入しても大丈夫.
すると
\[
p(\sigma)\begin{pmatrix}
a_n\\
b_n\\
c_n
\end{pmatrix}
=\tilde{A}(\sigma)(\sigma I-A)
\begin{pmatrix}
a_n\\
b_n\\
c_n
\end{pmatrix}
=\textbf{0}
\]となって,確かに作用素としての等式
\[
p(\sigma)=\det(\sigma I-A)=0
\]が得られる.なるほどね.
因みにa_n,b_nは役割が対称なので新たにd_n=a_n+b_nとおけば,
\[
\begin{pmatrix}
d_{n+1}\\
c_{n+1}
\end{pmatrix}
=
\begin{pmatrix}
3 & 4\\
2 & 3
\end{pmatrix}
\begin{pmatrix}
d_n\\
c_n
\end{pmatrix}
\]と書けるので,今度は
\[
\begin{vmatrix}
3-\sigma & 4\\
2 & 3-\sigma
\end{vmatrix}
=\sigma^2-6\sigma+1=0,
\]すなわち
\[
(\sigma^2-6\sigma+1)c_n=c_{n+2}-6c_{n+1}+c_n=0
\]なる漸化式が得られる.