LiNGAMモデル（因果推論）

因果推論周辺について調べた。

構造方程式モデルによる因果推論: 因果構造探索に関する最近の発展 http://www.ar.sanken.osaka-u.ac.jp/~sshimizu/papers/BSJ2012_Tutorial_final_web.pdf

非ガウス性を用いた線形非巡回なデータ生成過程部分の発見と同定 https://kaigi.org/jsai/webprogram/2012/pdf/258.pdf

ざっくり書くと、xを確率変数を縦にならべたベクトルとして,非巡回かつ潜在変数なしで、xを並び替えてBを下三角行列、eを非正規分布の確率ノイズとして x=Bx+e　 とかけるようなモデル。 結局、correlation-faithfulness の仮定がいるのか、チェックはどうするのか、eが非ガウスならOKなのかよく分からなかった。 とりあえず、Rのpcalgというパッケージの中のLINGAMで試してみた。

#1から2
x1=runif(N)-0.5
x2=x1*0.7+runif(N)-0.5
LINGAM(cbind(x1,x2))

$B
         [,1] [,2]
[1,] 0.000000    0
[2,] 1.582669    0

$Adj
      [,1]  [,2]
[1,] FALSE  TRUE
[2,] FALSE FALSE
#2から1
x2=runif(N)-0.5
x1=2*x2+runif(N)-0.5
LINGAM(cbind(x1,x2))

LINGAM(cbind(x1,x2))
$B
     [,1]     [,2]
[1,]    0 1.738995
[2,]    0 0.000000

$Adj
      [,1]  [,2]
[1,] FALSE FALSE
[2,]  TRUE FALSE

#1から2へガウスノイズ
x1=runif(N)-0.5
x2=x1*1.5+rnorm(N)
LINGAM(cbind(x1,x2))

$B
         [,1] [,2]
[1,] 0.000000    0
[2,] 1.845528    0

$Adj
      [,1]  [,2]
[1,] FALSE  TRUE
[2,] FALSE FALSE

ガウスでも多少は当たる！？ (上の例で160/200回くらいの正解率）

モデル妥当性のチェックも機械学習のcross-validationみたいなのが無いので、良く分からない。

• ブートストラップ法によるアプローチ: – サンプルサイズが小さいか非ガウス性が小さいかすれば、ブートスト ラップ標本についてのLiNGAMの結果は大きくばらつくはず

ふーむ。。。。

２声の自動生成

４声は悲惨だが２声の対旋律付けなら酷いレベルではない・・・とおもう。 ドフェミファソファミレの繰り返しのバスに対してソプラノの対旋律を生成した（両方生成する事もできる）。

www.youtube.com

G7でファ→ミの限定進行など和声っぽい音の動きが見られた。

ハ長調ではなくて旋法（ドリアとか）だと面白いのだけど、印刷物含めデータが多くないので、いずれ。

f:id:biones:20160307215903j:plain

上の画像の台形の音の音程について、(前の和音のMIDI距離,ソプラノの増分,テノールの増分)を特徴量としている。 この例だと(17,1,2)というベクトル。 探索方法はとりあえずで、焼きなまし法。 ドレミファソラシドという旋律とドの初期音からこのスケールを導くためにはひとつ前の和音も考慮すればいい（未実装）。

４声にしたときに、縦の響きについては流石に（台形の特徴量で）線形だとダメで（この特徴量だとG７の時は２度がOKみたいな状況からXORパターンがでてくると思われるから）、

(1)台形を一つふやす(4C3通り）
(2)カーネル法を使う

もっそり実験しようと思います。

3/23追記 1classSVM(と各種カーネル)でも4声は微妙。特徴量の次元が台形だけでも510次元,手持ちの課題を全部入力しても500和音程度にしかならないので（6000円の解答集を買えば3000和音くらいにはなる？）、実際の楽譜を使えるような方法に切り替えたほうがいいのだろうか。 転移音（和声外の音）や、4声（+旋律）→5声→3声みたいに特徴量の次元が不定で、そういうデータを上手く扱う方法はないのだろうか（ノンパラメトリックベイズとか、RNN？）。

kernelPCA(500→10〜200次元くらいまで落とす）→kmeans　でも微妙。教師なし学習の検討はできないので、相変わらず生成して私の耳で検証（楽で良いのだが）。異常検知の本とか買うと良いのだろーか。。。 いや、多分和声の解答集を買う事ですね。。。。

音楽データ分析(2)

島岡氏の和声の本のサンプルやら,図書館で借りてきた本の解答を手動でノーテーションソフトに入力して,music21等を使い、４声になっていない部分を補完。

和音の声部の横の流れはとりあえず考えないで、入力した音の基準で和音の縦の響きの「良さ」を計算してみる。 最大エントロピー法で、Xを和音,x_i∈Xを音程と(鍵盤の距離で0-34して
{p(X)=\frac{e^{\sum_{k} \omega * f_k(x) }}{Z} }
というモデルを考える。 素性fは音程があれば１そうでなければ０。重みωは３５次元で、L2ノルムの制約つき(ω2=1)。 和音に出てくる２音の音程（つまり４音なら4C2=6パターン）のXω/(２度のパターン数)をスコアとしました。 入力した和音は4音350くらい。全部C-major。

f:id:biones:20160219174009p:plain 音程のスコアの分布。x軸はMIDIの距離(オクターブ→12)

無事オクターブ(12),５度(7),3度(3,4)あたりが高くなった。

本当は音程ではなく音そのものでやりたかったのだけど(低いド-ミと高いド-ミの意味はだいぶ違うと思うから）、思っていたよりかなりデータが手に入りにくいので、低次元な方向でやりたいと思う。 次は隣り合った２和音について、「台形」の(前の和音の２声部の縦の音程、下の声部が変化分,上の声部の変化分)を素性としてやりたいと思う。これだと連続・平達５度等が考慮される・・・ハズ。

XTを教師データの音,x0∈XTを探索を開始する和音,Xを探索パラメータωに基づいて生成された音, scoreωをXに対するωで評価したスコアとして scoreω(X(x0,ω))とscore_ω(XT)の差をlossとしてωを最適化する、探索を組み込んだ将棋のボナンザメソッドみたいな事できないだろーか。。。

転調は調を隠れ変数、転移音は転移音用の素性を用意すれば、大分実際の曲のデータを使えるようになる・・・と思う。

pandas,numpy.matrix,listの行き来はまだまだ慣れないですorz

音楽データ分析

昔、マルコフ連鎖インベンションみたいな事をやったり、ルールベースの対位法もどきみたいな事をやって放置していたけどようやく再開。

やり方は,楽譜ソフトに音を入力して music21: a Toolkit for Computer-Aided Musicology というPythonのライブラリで読み込み、足りない音の補完など形を整えてcsvで出力→Rで解析。

ちなみにmusic21の日本語の紹介記事 [楽譜][可視化]Pythonで音楽学 d.hatena.ne.jp

目的は、4声帯（大体）の音の予測。 i番目の和音を{C_i}として、{C_i-1,C_i-2,C_i-{y_i}}から目的の音（ソプラノとか）{y_i}を予測するモデル。

線形分離とは無縁だと思われるのでとりあえずてきとーにランダムフォレストしてみたらC-majorのみの５０和音程度で(島岡和声の１巻より)error rate27%くらい。これはデータを増やせば行けるのでは！？

隠れマルコフモデルで、コード（とか緊張度とか）を潜在変数として、具体的な音はその実現値とする思想がしっくりきていたけど、クラシックだと何となく違うような気もするのですお。iステップ目のサイコロの出目からi+1ステップ目の和音選択への影響が大きい、とかむしろ調もコードも無しでいいのかなぁ・・・とか。

例えば　　 ω: C→Dm7→G7
tk1→t(k+1)1→t(k+2)1・・・
tk2
tk3
tk4
↑ kステップ目の音,声部(max4)

みたいなモデルで
P(t_k1 | t_k-1,H_k-1,ω_k-1)　 H_k={ ∪i t_ki}

P(ωk|ωk-1)　　コードの遷移
P(t_k2|t_k1,H_k-1,ω_k-1) 　音の出力
みたいな。

自己回帰とか多段のマルコフモデルもある様ですね。

ベートーヴェンとスクリャービンの楽譜を混ぜて生成した音を使って作った曲をドヤ顔で弾くのが夢。

2種類の季節変動をstanで推定

状態遷移の式 α(t+1)=T*αt+ε　　(記号はKFASから） でTの中に2種類のQが入るだけなので、何も問題がないはずです。

nn=3
a=rep(c(1,2,3),5*nn)
b=rep(c(5,5,0,0,0),3*nn)
d=10+a+b
ts.plot(d[1:30])

f:id:biones:20151019181030j:plain

scode="
data{
  int N;
  real y[N];
  
}
parameters{
  real sigma;
  real trend[N];
  real s_trend;

  real season[N];
  real s_season;

  real season2[N];  #2つ目の季節変動
  real s_season2;
}

model{
  real ssum;
  real mu[N];
  
  for(t in 2:N){
    trend[t]~normal(trend[t-1],s_trend);
  }

  for(t in 5:N){
    ssum<-0;
    for(k in 1:4){
      ssum<-ssum+season[t-k];
    }
    season[t]~normal(-ssum,s_season);
  }

  for(t in 3:N){
    ssum<-0;
    for(k in 1:2){
      ssum<-ssum+season2[t-k];
    }
    season2[t]~normal(-ssum,s_season2);
  }

  for(t in 1:N){
    mu[t]<-season[t]+trend[t]+season2[t];
    y[t]~normal(mu[t],sigma);
  }

}

"

library(rstan)
N=length(d)
slist=list(y=d,N=N)
fit=stan(model_code = scode,data=slist ,iter=1000,warmup=300,chains=1)

#stanの後処理
res=stan.ss(fit,y=d,varnames = c("season","season2"))

r=d-(res$trend+res$season+res$season2)
#> mean(r^2)
#[1] 3.102593e-21

f:id:biones:20151019180832j:plain (yと推定結果が重なってます)

f:id:biones:20151019180812j:plain