Julia #4: Nonnegative Matrix Factorization

————————————
ベイズ推論による機械学習入門　須山敦志　著　を読む
————————————
NMF：nonnegative matrix factorization
————————————
まずは、Anacondaでscipyをインストールしておこう。
organ.wavという学習用データファイルを使用するので、ここではパスを間違えないように絶対パスで記載する。

まずは、NMFモジュールとして、NMF.jlのコードを一気にJupyterへドンと実行。

"""
Variational inference for Bayesian NMF
"""
module NMF
using Distributions
using StatsFuns.logsumexp, SpecialFunctions.digamma

export NMFModel
export sample_data, VI

途中省略.....

function VI(X::Array{Int64, 2}, model::NMFModel, max_iter::Int)
    K = size(model.a_t, 2)
    D, N = size(X)
    S, A_t, B_t, A_v, B_v = init(X, model)
    for iter in 1 : max_iter
        # latent
        S = update_S(X, A_t, B_t, A_v, B_v)
        A_v, B_v = update_V(S, A_t, B_t, model)

        # param
        A_t, B_t = update_T(S, A_v, B_v, model)
    end
    
    return update_model(A_t, B_t, model), S, A_t.*B_t, A_v.*B_v
end

end

"""

Variational inference for Bayesian NMF

"""

module NMF

using Distributions

using StatsFuns.logsumexp, SpecialFunctions.digamma

export NMFModel

export sample_data, VI

途中省略.....

function VI(X::Array{Int64, 2}, model::NMFModel, max_iter::Int)

K = size(model.a_t, 2)

D, N = size(X)

S, A_t, B_t, A_v, B_v = init(X, model)

for iter in 1 : max_iter

# latent

S = update_S(X, A_t, B_t, A_v, B_v)

A_v, B_v = update_V(S, A_t, B_t, model)

# param

A_t, B_t = update_T(S, A_v, B_v, model)

end

return update_model(A_t, B_t, model), S, A_t.*B_t, A_v.*B_v

end

Out: NMF

using NMF

using NMF

次にdemo_NMF.jlを、動作を確認しながら、少しずつ入力実行していく。その歳、moduleのソースコードや、データファイルのパスには要注意（ここでは絶対パスにした）。

##############################
## Audio decomposition demo using NMF

using PyPlot, PyCall
using DataFrames
using Distributions

push!(LOAD_PATH, "/Users/******/Desktop/BayesBook-master/src")
import NMF
@pyimport scipy.io.wavfile as wf

##############################

## Audio decomposition demo using NMF

using PyPlot, PyCall

using DataFrames

using Distributions

push!(LOAD_PATH, "/Users/******/Desktop/BayesBook-master/src")

import NMF

@pyimport scipy.io.wavfile as wf

# load data
wavfile = "/Users/*******/Desktop/BayesBook-master/data/organ.wav"
fs, data = wf.read(wavfile)

# load data

wavfile = "/Users/*******/Desktop/BayesBook-master/data/organ.wav"

fs, data = wf.read(wavfile)

Out: (44100, Int16[0 0; 0 0; … ; 0 0; 0 0])
WAVファイルの情報はPyObject Arrayであるdataに取り込まれている、
scipyのspecgram()でフーリエ変換。ここのdata[10000:318000,2]が、Julia ver 0.7ではにっちもさっちもブルドッグ。ver.0.6.4に下げて、直接、Jupyterで入力すれば何のことはない。

figure("data")
clf()
Pxx, freqs, t, pl = specgram(data[10000:318000,2], Fs=fs, NFFT=256, noverlap=0)
xlabel("time [sec]")
ylabel("frequency [Hz]")
ylim([0,22000])

figure("data")

clf()

Pxx, freqs, t, pl = specgram(data[10000:318000,2], Fs=fs, NFFT=256, noverlap=0)

xlabel("time [sec]")

ylabel("frequency [Hz]")

ylim([0,22000])

Out: (0, 22000)

# model
D, N = size(Pxx)
K = 2
a_t = 1.0
b_t = 1.0
a_v = 1.0
b_v = 100.0
prior = NMF.NMFModel(a_t*ones(D,K), b_t*ones(D, K), a_v, b_v)

# model

D, N = size(Pxx)

K = 2

a_t = 1.0

b_t = 1.0

a_v = 1.0

b_v = 100.0

prior = NMF.NMFModel(a_t*ones(D,K), b_t*ones(D, K), a_v, b_v)

Out: NMF.NMFModel([1.0 1.0; 1.0 1.0; … ; 1.0 1.0; 1.0 1.0], [1.0 1.0; 1.0 1.0; … ; 1.0 1.0; 1.0 1.0], 1.0, 100.0)

# inference
max_iter = 100
posterior, S_est, T_est, V_est = NMF.VI(Int64.(round.(Pxx)), prior, max_iter)
X = T_est * V_est

# inference

max_iter = 100

posterior, S_est, T_est, V_est = NMF.VI(Int64.(round.(Pxx)), prior, max_iter)

X = T_est * V_est

Out:129×1203 Array{Float64,2}:
0.758165 0.345822 1.55558 … 0.547368 0.315923 1.79641
4.70464 2.14591 9.65305 3.3964 1.95962 11.1472
1.88055 0.857836 3.85757 1.35844 0.786708 4.4563
0.119958 0.062215 0.132823 0.181609 0.444062 0.344534
0.0132027 0.00667184 0.0172718 0.0177633 0.0396456 0.0365076
0.00425499 0.00212669 0.00592189 … 0.0054267 0.0115406 0.0115765
0.0341117 0.0180025 0.0330729 0.0555813 0.142612 0.100473
1.39303 0.738848 1.2951 2.31633 6.01692 4.13258
2.1343 1.13212 1.98252 3.55037 9.22472 6.33255
1.16528 0.617831 1.08671 1.93482 5.02155 3.45515
0.38544 0.204349 0.359622 … 0.639837 1.66039 1.14277
0.563663 0.299005 0.523379 0.937812 2.43693 1.67252
0.520796 0.276265 0.483573 0.866491 2.2516 1.54532
? ? ?
1.01537e-6 5.35847e-7 9.84728e-7 1.65421e-6 4.24402e-6 2.99053e-6
1.01537e-6 5.35847e-7 9.84728e-7 1.65421e-6 4.24402e-6 2.99053e-6
1.01537e-6 5.35847e-7 9.84728e-7 1.65421e-6 4.24402e-6 2.99053e-6
1.01537e-6 5.35847e-7 9.84728e-7 … 1.65421e-6 4.24402e-6 2.99053e-6
1.01537e-6 5.35847e-7 9.84728e-7 1.65421e-6 4.24402e-6 2.99053e-6
1.01537e-6 5.35847e-7 9.84728e-7 1.65421e-6 4.24402e-6 2.99053e-6
1.01537e-6 5.35847e-7 9.84728e-7 1.65421e-6 4.24402e-6 2.99053e-6
1.01537e-6 5.35847e-7 9.84728e-7 1.65421e-6 4.24402e-6 2.99053e-6
1.01537e-6 5.35847e-7 9.84728e-7 … 1.65421e-6 4.24402e-6 2.99053e-6
1.01537e-6 5.35847e-7 9.84728e-7 1.65421e-6 4.24402e-6 2.99053e-6
1.01537e-6 5.35847e-7 9.84728e-7 1.65421e-6 4.24402e-6 2.99053e-6
1.01537e-6 5.35847e-7 9.84728e-7 1.65421e-6 4.24402e-6 2.99053e-6

# visualize
figure("T")
clf()
for k in 1 : K
    subplot(K,1,k)
    plot(T_est[:,k], linewidth=1.0)
    xlim([0, D])
    ylim([0, ylim()[2]])
end

# visualize

figure("T")

clf()

for k in 1 : K

subplot(K,1,k)

plot(T_est[:,k], linewidth=1.0)

xlim([0, D])

ylim([0, ylim()[2]])

end

figure("V")
clf()
for k in 1 : K
    subplot(K,1,k)
    plot(V_est[k,:], linewidth=1.0)
    xlim([0,N])
    ylim([0, ylim()[2]])
end
show()

figure("V")

clf()

for k in 1 : K

subplot(K,1,k)

plot(V_est[k,:], linewidth=1.0)

xlim([0,N])

ylim([0, ylim()[2]])

end

show()

Science To Medicine

Just My Daily Study Note by ts.anesth.kpum

Julia #4: Nonnegative Matrix Factorization