2声の自動生成

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

www.youtube.com

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

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

f:id:biones:20160307215903j:plain

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

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

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

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

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

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

音楽データ分析(2)

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

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

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

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

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

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のみの50和音程度で(島岡和声の1巻より)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

UCI Machine Learning RepositoryよりPoker Hand

UCI Machine Learning Repository: Poker Hand Data Set

ポーカーの手の判定。プログラミングの課題でよくあるけど、統計的な分類で判定します。

> head(df)
  C 1 C 2 C 3 C 4 C 5 C 6 C 7 C 8 C 9 C 10 C 11
1   3  12   3  11   3  13   3  10   3    1    9
2   4  10   4  11   4   1   4  13   4   12    9
3   4   1   4  13   4  12   4  11   4   10    9
4   1   2   1   4   1   5   1   3   1    6    8
5   1   9   1  12   1  10   1  11   1   13    8
6   2   1   2   2   2   3   2   4   2    5    8

0: Nothing in hand; not a recognized poker hand 
1: One pair; one pair of equal ranks within five cards 
2: Two pairs; two pairs of equal ranks within five cards 
3: Three of a kind; three equal ranks within five cards 
4: Straight; five cards, sequentially ranked with no gaps 
5: Flush; five cards with the same suit 
6: Full house; pair + different rank three of a kind 
7: Four of a kind; four equal ranks within five cards 
8: Straight flush; straight + flush 
9: Royal flush; {Ace, King, Queen, Jack, Ten} + flush 


役の分布 N=25008 
0     1     2     3     4     5     6     7     8     9 
12493 10599  1206   513    93    54    36     6     5     3 

ロイヤルフラッシュは0.0001199616 !1/一万

元の特徴量はカードのスートと数で、1-9の役を目的変数として(C11,当然質的変数に直して)、そのままやっても上手くいくハズはないので、 1.カードの重複枚数の1番大きい数 2.カードの重複枚数の2番目に大きい数 3.スートが5つ揃っていらTrueそうでなければFalse 4.ストレートの判定 5.ストレートフラッシュの判定

と特徴量を直して、randomForestで分類。

こんなに追加したら普通に判定するのと大して変わらないような気もするのですが、該当部分は23行で書けているのでいいのかなと。モノ的にエラー率0%じゃないとダメでしょうし。 forの中でソートを2回しているのはスルーしてくださいw

df=read.csv("https://archive.ics.uci.edu/ml/machine-learning-databases/poker/poker-hand-training-true.data")
colnames(df)=paste("C",as.character(1:11))


#特徴量の追加
ncard1=c()
ncard2=c()
suit=c()
straight=c()
sf=c()
for(i in 1:N){
  dsuit=df[i,c(1,3,5,7,9)]
  d=df[i,c(2,4,6,8,10)]
  hoge=sort(table(as.numeric(d)),decreasing = T)
  d1=hoge[1]
  d2=hoge[2]
  ncard1=c(ncard1,d1)
  ncard2=c(ncard2,d2)
  #スート
  if(length(unique(as.numeric(dsuit)))==1){
    same_suit=1
  }else{
    same_suit=0
  }
  suit=c(suit,same_suit)
  
  d=sort(d)
  if(d[5]-d[1]==4 || setequal(d,c(1,10,11,12,13))){
    st=T
  }else{
    st=F
  }
  straight=c(straight,st)
  hoge=ifelse(max(d)==13 && min(d)==1 && st,T,F)
  sf=c(sf,hoge)
}

df2=cbind(ncard1,ncard2,suit,straight,sf)
df2=data.frame(df2)
df2$straight=factor(df2$straight)
df2$sf=factor(df2$sf)

df2=cbind(df2,y=df[,11])

df2$y=factor(df2$y)
df2$suit=factor(df2$suit)
rownames(df2)=as.character(1:N)


N=nrow(df2)
trn=sample(N,floor(N)*0.7)
train=df2[trn,]
test=df2[-trn,]


library(randomForest)
tuneRF(y = train[,6],x=train[,-6],type="classification")

rf=randomForest(y~.,data=train,mtry=5,type="classification")

pred=predict(rf,newdata=test)
cm=table(test$y,pred)
print(cm)

   pred
       0    1    2    3    4    5    6    7    8    9
  0 3824    0    0    0    0    0    0    0    0    0
  1    0 3121    0    0    0    0    0    0    0    0
  2    0    0  352    0    0    0    0    0    0    0
  3    0    0    0  156    0    0    0    0    0    0
  4    0    0    0    0   18    0    0    0    0    0
  5    0    0    0    0    0   17    0    0    0    0
  6    0    0    0    0    0    0    8    0    0    0
  7    0    0    0    0    0    0    0    4    0    0
  8    0    0    0    0    0    0    0    0    2    0
  9    0    0    0    0    0    0    0    0    0    1

グリッドの人の移動と軌跡のクラスタリング

Pythonjavascript

f:id:biones:20150922194909p:plain

http://biones.dip.jp/space/space.py (鯖が止まってると動かない可能性あり)

的なモノを昔作って何か分析しようとおもいつつ放置していたのに再着手。 青色の点が目的地。

歩行軌跡間の「近さ」はデータの長さが違っても類似度が計算できるDynamic Time Warping距離

ダイナミックタイムワーピング距離に基づくストリーム処理 http://www.ieice.org/iss/de/DEWS/DEWS2007/pdf/l6-5.pdf

で測って、距離行列を作成、k-medroid(pamパッケージ)でクラスタリングして、何処に向かっているか判定します。

生成の方のコードは長くかつ汚いので省略(このRコードもきたな・・・)。

m=read.csv("~/Desktop/traject.csv")
m=m[,-1]

library(dtw)

hns=sort(m$human_no)
hns=unique(hns)

library(ggplot2)
gdf=data.frame(matrix(ncol = 4)) #ggplot用
colnames(gdf)=c("hn","x","y","cl")
#gdf=gdf[-1,]

#プロット用
plt=function(gdf){
  gdf$cl=factor(gdf$cl)
  for(hn in hns){
    x=subset(m,human_no==hn)$x
    y=subset(m,human_no==hn)$y
    
    hoge=pm$clustering
    cl=hoge[names(hoge)==as.integer(hn)]
    if(length(cl)<=0){
      next
    }
    tdf=cbind(hn,x,y,cl)
    gdf=rbind(gdf,tdf)
  }
  gp=ggplot(data=gdf,aes(x=x,y=y,color=cl))+geom_point()
  
  print(gp)
}


#距離行列の計算

hoge=combn(hns,2)
dmat=data.frame(NULL)
cc=combn(hns,2)
for(j in 1:dim(hoge)[2]){
  c=cc[,j]
  x1=subset(m,human_no==c[1])$x
  y1=subset(m,human_no==c[1])$y
  
  x2=subset(m,human_no==c[2])$x
  y2=subset(m,human_no==c[2])$y
  
  dist=dtw(x=cbind(x1,y1),y=cbind(x2,y2))  #DTW距離
  
  dmat[paste(c[1]),paste(c[2])]<-dist$distance
  if(j%%1000==0){
    print(j) 
  }
}

for(i in 1:dim(dmat)[1]){
  for(j in 1:dim(dmat)[2]){
    if(is.na(dmat[i,j])){
      dmat[i,j]<-0
    }
  }
}

#クラスタリング
pm=pam(dmat,k=5)#クラスタ数は既知
pred=pm$clustering

mokutekiti=floor(as.numeric(names(pred))/1000)+1 #human_noの4桁目で目的地を分けているので
res=table(mokutekiti,pred)
print(res)

plt(gdf)

目的地の生成部分

f.generate("S1",10,10)
f.generate("S2",70,60)
f.generate("S3",10,40)
f.generate("S4",10,20)
f.generate("S5",80,70) 

結果。

          pred
mokutekiti   1   2   3   4   5
         1 100   0   0   0   0
         2   0 100   0   0   0
         3   0   3  94   3   0
         4   0   0  17  83   0
         5   0   0   0   4  95
  

f:id:biones:20150922202413j:plain

ggplot頑張ったけど見づらいですorz

距離行列の計算が重くて500人程度でも時間がかかるので、隠れマルコフモデルで再チャレンジしてみようかと・・・。