No science, No life.

もっと科学を身近に

カオス

「部屋がカオス」といえば部屋が乱雑に散らかっていることを指すと思いますが、「カオス理論」というときのカオスは単に乱雑ということとはもちろん違います。

有名なロジスティック写像を例にしてカオスの様子を覗いてみましょう。ロジスティック写像は、

{ x_{n+1} = a x_n(1-x_n)}

であらわされる単純な漸化式です。電卓でも簡単に計算できますね。たとえば、定数a=0.5, 初期値x(1)=0.9とすれば、x(2) = 0.5*0.9*(1-0.9) = 0.045, x(3) = 0.5*0.045*(1-0.045) = 0.0215, ...という具合です。

一見したところ単純で、しかも決定論的 (前の値によって次の値が完全に決まる) なこの式ですが、定数aの値によって多彩な振る舞いをします。aの値を変えて様子を見てみましょう。

下のMatlabスクリプトで、まず、a=0.5のときの様子を描いてみます。初期値x(1)は適当に0.9としました。

clear;
a = 0.5;
x(1) = 0.9;

for n=1:99
    x(n+1) = a * x(n) * (1-x(n));
end

figure(1); clf
plot(1:100,x,'.-')
xlabel('n')
ylabel('x(n)')
title(strcat('a=',num2str(a)))
ylim([0 1])

f:id:neocortex:20161120170602p:plain

a=0.5だと値がすぐに0に収束します。では次にa=2として描いてみます。

f:id:neocortex:20161120171004p:plain

この場合も特定の値に収束しました。次にa=3.3としてみましょう。

f:id:neocortex:20161120171238p:plain

この場合は収束せず、2つの値のあいだで振動しています。では、a=3.5はどうでしょうか。

f:id:neocortex:20161120172309p:plain

この場合も振動していますが、値が4通りあることがわかります。では、a=3.8としてみます。

f:id:neocortex:20161120195011p:plain

すると先ほどまでとは違い、さまざまな値をとりながら (一見したところ) 不規則に揺れ動いています。このような状態がカオスです。さらにa=4とすると、0から1のすべての値をとるカオスとなります。

f:id:neocortex:20161120195846p:plain

以上のような状態の変化を見通しよくするために、下のスクリプトで「分岐図」を描いてみます。aの値を0から4のあいだで小刻みに変えて、x(100)まで計算し、最後の20点のみをプロットしています (例として挙げたaの値には赤線を引きました)。

clear;
x(1) = 0.9;
figure(1); clf 
for a = 0:0.001:4
    for ii=1:99
        x(ii+1) = a * x(ii) * (1-x(ii));
    end
    plot(a*ones(1,20),x(81:100),'b.');
    hold on
end
title('分岐図')
xlabel('a')
ylabel('x(n)')

plot([0.5 0.5],[0 1],'r:')
plot([2.0 2.0],[0 1],'r:')
plot([3.3 3.3],[0 1],'r:')
plot([3.5 3.5],[0 1],'r:')
plot([3.8 3.8],[0 1],'r:')
plot([4.0 4.0],[0 1],'r:')

f:id:neocortex:20161120201357p:plain

f:id:neocortex:20161120205752p:plain

上のプロットを見ると、aの増加にともなって、値の収束が1点, 2点, 4点...と増えていることが分かります。これを「周期倍分岐」といいます。そして、a=3.56付近以上でカオスが登場することが分かります。

分岐図をよく見てみると、同じような構造が何回も登場することが分かります。下のスクリプトでa=3.8から3.9の範囲をこまかくプロットしなおしてみましょう。

clear;
x(1) = 0.9;
figure(1); clf 
for a = 3.8:0.0001:3.9
    for ii=1:199
        x(ii+1) = a * x(ii) * (1-x(ii));
    end
    plot(a*ones(1,20),x(181:200),'b.');
    hold on
end
title('分岐図')
xlabel('a')
ylabel('x(n)')
xlim([3.8 3.9])

f:id:neocortex:20161120211341p:plain

f:id:neocortex:20161120211408p:plain

a=3.84から3.85のあたりを拡大してみるとここにも分岐が見えます。拡大してもおなじ形が見えることを「自己相似」といい、完全な自己相似性をもつ図形は「フラクタル」と呼ばれます。

さて、カオスの重要な特徴として、将来の予測が (事実上) 不可能だという点があります。現実の世界の数値にはつねに誤差があります (たとえば体重を測るにしても測定誤差があり、無限の精度で測定するということはできません)。カオスな世界では初期値のわずかな誤差があっという間に増幅されてしまい、遠い将来の値を予測することが事実上できません。これを「初期値鋭敏性」と呼びます。

下のスクリプトで初期値鋭敏性をシミュレーションしてみましょう。カオスな条件 (a=4) で、初期値が0.9の数値列x(n)と、初期値が0.9から1/10万だけずれた数値列y(n)を比較してみます。

% ロジスティック写像の計算
clear;
a = 4;
x(1) = 0.9;
y(1) = x(1)*1.00001;

for n=1:99
    x(n+1) = a * x(n) * (1-x(n));
    y(n+1) = a * y(n) * (1-y(n));
end

% x(n), y(n) をプロット
figure(1); clf
subplot(3,1,1:2)
plot(1:100,x,'b.-',1:100,y,'r.-')
ylabel('x(n), y(n)')
title(strcat('a=',num2str(a)))
ylim([0 1])

% x(n)-y(n) をプロット
subplot(3,1,3)
plot(1:100,x-y,'k.-')
xlabel('n')
ylabel('x(n)-y(n)')

f:id:neocortex:20161120215043p:plain

n=14あたりまでは、青線のx(n)と赤線のy(n)はほとんど同じ値でプロットが重なっていますが、それ以降はまったく別の値をとり、プロットが重ならないことが分かります。1/10万の誤差がいつの間にか増幅されたわけです。

以上のロジスティック写像にはランダムなノイズはまったく含まれておらず、完全に決定論的です。このように、「決定論的なのに将来が予測できない」というところにカオスの衝撃があります。

カオス力学の基礎

カオス力学の基礎

「吾輩は猫である」のテキストを分析してみる (4)

「吾輩は猫である」のテキストを分析してみる (1)
「吾輩は猫である」のテキストを分析してみる (2)
「吾輩は猫である」のテキストを分析してみる (3)


前回、何も工夫せずに頻出単語を取りだすと、記号や助詞・助動詞ばかりがでてきました。作品の特徴を捉えた頻出単語帳を作ろうと思うと、品詞を絞ったほうがよさそうです。そこで今回は、名詞 (の一部) のみの頻出単語を並べてみました。

f:id:neocortex:20160910230659p:plain

吾輩・寒月・先生・猫・細君などの単語が見えます。「吾輩は猫である」の特徴を捉えた単語帳になってきたでしょうか。コードは下の通りです。

% 形態素分析結果 result の読込
load('neko.mat')

% 名詞を抽出
idx1 = strcmp(result(:,2),'名詞');
idx2 = strcmp(result(:,3),'非自立');
idx3 = strcmp(result(:,3),'数');
result_noun = result(idx1 & ~idx2 & ~idx3,:);

% セル配列からカテゴリカル配列を作成
cat_array = categorical(result_noun(:,1));

% 単語カテゴリーを作成
category = categories(cat_array);

% 各カテゴリの出現回数
counts = countcats(cat_array);
% 出現回数を出現パーセントに変換
prct = counts/sum(counts)*100;

% 出現パーセントに従って降順でソート
[~,I] = sort(prct,'descend');
prct_sorted = prct(I);
category_sorted = category(I);

% グラフ表示
figure(1); clf
x = 1:length(prct_sorted);
y = prct_sorted;
plot(x,y,'o-')
xmax = 20;
axis([0 xmax 0 prct_sorted(1)+0.2])
text(x(1:xmax),y(1:xmax)+0.1,category_sorted(1:xmax))
title('「吾輩は猫である」形態素分析')
xlabel('名詞出現ランク')
ylabel('名詞出現パーセント')

吾輩は猫である (新潮文庫)

吾輩は猫である (新潮文庫)

「吾輩は猫である」のテキストを分析してみる (3)

「吾輩は猫である」のテキストを分析してみる (1)
「吾輩は猫である」のテキストを分析してみる (2)



前回までで、「吾輩は猫である」を形態素分析にかけたセル配列が手に入りました。今回は、単語の出現頻度を調べてみようと思います。前回のコードで、

in = txt(1:1000);

として、最初の1000文字だけを取り出し、下のコードで出現頻度をプロットしてみました。

% セル配列からカテゴリカル配列を作成
cat_array = categorical(result(:,1));

% 単語カテゴリーを作成
category = categories(cat_array);

% 各カテゴリの出現回数
counts = countcats(cat_array);
% 出現回数を出現パーセントに変換
prct = counts/sum(counts)*100;

% 出現パーセントに従って降順でソート
[~,I] = sort(prct,'descend');
prct_sorted = prct(I);
category_sorted = category(I);

% データを保存
save('neko.mat','result','category_sorted','prct_sorted')

% グラフ表示
figure(1); clf
x = 1:length(prct_sorted);
y = prct_sorted;
plot(x,y)
title('「吾輩は猫である」形態素分析')
xlabel('出現ランク')
ylabel('出現パーセント')
xmax = max(x);
axis([0 xmax 0 prct_sorted(1)+1])

結果は下のようになりました。もとの単語数が少ないのでグラフがガタガタしています。

f:id:neocortex:20160828225741p:plain

次に、全文を突っ込んでみると次のようになりました。軸に張りついていて分かりにくいですが、単語数が増えてグラフが滑らかになりました。

f:id:neocortex:20160828230123p:plain

同じデータを、対数グラフでプロットしてみると、下のように広い範囲でほぼ直線になりました。

loglog(x,y)

f:id:neocortex:20160828230623p:plain

対数グラフ上で直線になる関係を冪乗則 (べきじょうそく) といいます。これは、単語の出現頻度分布をふくめて、自然界のさまざまな状況で現れる関係性だそうです。「吾輩は猫である」の文章も例外ではなく、冪乗則に従うことがわかります。

さて、上位に出現するのはどんな言葉でしょうか?下のようにトップ50位を表示してみました。

plot(x,y,'o-')
xmax = 50;
axis([0 xmax 0 prct_sorted(1)+1])
text(x(1:xmax),y(1:xmax)+0.5,category_sorted(1:xmax))

f:id:neocortex:20160828232259p:plain

すると、当然と言えば当然ですが、最上位は「の」「。」「て」「、」「は」…となっていて、文章に必ず登場するような助詞・助動詞や記号ばかりです。意味のある名詞としては「吾輩」がかなり多いはずと思い見てみたら、47位にランクインしていました。

xlim([45 50])

f:id:neocortex:20160828232829p:plain

文章を特徴付ける単語を抽出してこようと思うと、記号・助詞・助動詞などは除いて名詞などに焦点を当てたほうが良さそうです。次回以降、やってみようと思います。

吾輩は猫である (新潮文庫)

吾輩は猫である (新潮文庫)

「吾輩は猫である」のテキストを分析してみる (2)

前回のneko.txtをMeCabの分析にかけていきます。

コードはこちら。neko.txtの全文をいきなり投入すると時間がかかるので、ここでは最初の100文字だけにしています。

% 「吾輩は猫である」のテキストを読込
fileID = fopen('neko.txt');
txt = fread(fileID,'*char')';
fclose(fileID);

% MeCab DLLを読込
[notfound,warnings] = loadlibrary('libmecab.dll','mecab.h');
mecab = calllib('libmecab','mecab_new',1,libpointer('stringPtrPtr',{'MeCab'}));

% MeCabで形態素分析
in = txt(1:100);
mecab_result = calllib('libmecab','mecab_sparse_tostr',mecab,in);

そうすると下のような結果が出てきます。

f:id:neocortex:20160828010041p:plain

この分析結果は、ズラズラッとつらなったテキスト形式でこのままでは扱いにくいので、下のようにセル配列に代入することにします。

ここで、MeCabの結果のフォーマットは「表層形\t品詞,品詞細分類1,品詞細分類2,品詞細分類3,活用型,活用形,原形,読み,発音」という10要素なのですが、改行記号などの場合は、strsplitで要素を分割したときに10要素に満たないことがあって扱いにくく、まぁ記号はいらないだろうということで消しています。

% MeCabの結果をセル配列に代入
row = strsplit(mecab_result,'\n');  % MeCabの結果をまず行ごとに分割
n = length(row)-2;                  % 行数を取得 (最後の2行は'EOS'等なので不要)
result = cell(n,10);                % 結果を格納する空配列を用意
idx = false(n,1);                   % あとで消去する行のインデックス
for ii=1:n
    % 空白と「,」で要素ごとに分割
    str = strsplit(row{ii},{'	',','});    
    if length(str)<10
        % 記号などの場合は要素数が10に満たないので消去
        idx(ii) = true;
    else
        % 結果を格納
        result(ii,:) = str;
    end
end
% いらない行を消去
result(idx,:)=[];

% DLLの開放
clear mecab
unloadlibrary('libmecab')

これで、下のような結果が得られました。

f:id:neocortex:20160828010928p:plain

「吾輩は猫である」のテキストを分析してみる (1)

青空文庫から夏目漱石の「吾輩は猫である」の全文をダウンロードして遊んでみます。テキストファイルをダウンロードして開いてみると下のようになっていました。

f:id:neocortex:20160827234426p:plain

ヘッダと (上の画像では見えていませんが) フッタと、本文中にたくさんルビなどの注釈がついています。これらはテキスト分析にはいらないので除きます。

ヘッダ・フッタは1箇所だけなのでメモ帳で手動で消しました。本文中の注釈は膨大な量なので、Matlab正規表現を使って消すことにしました。コードはこちら。

% load text data
fileID = fopen('wagahaiwa_nekodearu.txt');
txt = fread(fileID,'*char')';
fclose(fileID);

% remove annotations
pat = '《[^》]*》'; 
txt = regexprep(txt, pat, '');
pat = '[#[^]]*]'; 
txt = regexprep(txt, pat, '');
pat = '|'; 
txt = regexprep(txt, pat, '');

% save text data
fileID = fopen('neko.txt','w');
fprintf(fileID,txt);
fclose(fileID);

これで、下のように注釈がとれたテキストneko.txtができました。

f:id:neocortex:20160827235233p:plain

著作権の消滅した書籍を自由に読める「青空文庫」

著作権の消滅した書籍などを自由に読める「青空文庫」というサイトがあります。この記事を書いている時点で13,752の作品が登録されているそうで、誰もが知っているような文豪の作品がたくさんあります。

いろんな作品を無料で読めるという点もすごいですが、整った長文テキストデータという点でも魅力的です。ここからデータを拝借してすこし遊んでみようと思います。

形態素分析エンジンMeCabをMatlab (64bit)で使う

形態素分析エンジンMeCabMatlabから使ってテキストをササッと解析できるようにしたい!、ということでやってみました。ちょっと手間だったので書き残しておきます。

下のページにまさにこれを実現するためのツールがあるのですが、「現在配布されている MeCab の dll は 32bit version であるため、64bit 版のMATLABと一緒に使うためには、dll のビルドが必要となります」という不穏な忠告が。自分の場合は64 bit版Matlabなので、まさにこれにあたります…。

Simple Text Miner for Japanese - File Exchange - MATLAB Central

仕方がないので、dllをビルドするという自分には難易度の高い旅に出ました。

blogs.mathworks.com

上のブログ記事の中盤以降にある "What about Japanese Text?" という項目を探して、そのとおりにおこないます。

  1. Matlab 64 bit版を持っている
  2. Matlab に適合した (たぶんC++) コンパイラを持っている → 自分の場合は、無料のVisual studio 2013 Communityがインストールしてあったので、これを使いました。
  3. Takuya's instructionに従いなさい

という3ステップで、あら簡単、と思いきやTakuya's instructionがヘビーでした。以下、インストラクションのPDFを見つつ読んでください。

  1. MeCab 0.966のバイナリパッケージをインストール → インストール時に、文字のエンコーディングを聞いてきますが、デフォルトのSHIFT-JISで自分はうまくいきました。
  2. MeCab 0.966のソースファイルをダウンロード&解凍 → Windows環境では.tar.gzファイルを解凍するにはソフトウェアが必要です。自分はLhaplusを使いました。
  3. 解凍されたなかの"src"というフォルダを見つけます。
  4. このフォルダのなかにあるいろんなファイルの内容を、インストラクションの通りに修正します。
  5. Visual studioコマンドプロンプト (いろいろあるが、自分はVS2013 x64 Native Tools Command Promptでうまくいきました) を起動します。
  6. フォルダをさきほどの"src"にもっていきます。フォルダを変更するコマンドはcd。
  7. "Make"とタイプ。
  8. すると、ずらずらっとビルドが始まり、"libmecab.dll"が生成されます。
  9. このdllファイルと、"mecab.h"ファイルを、Matlabスクリプトを置くフォルダにコピー。これで、MatlabMecabを使う準備完了です。ふぅ。

では、Matlabを起動してサンプルコードをコピペしてみます。まずは下のコードでDLLを読み込み。僕の環境では"警告: データ型 'mecab_node_tPtr' (構造体 mecab_node_t で使用) は存在しません…"というエラーが出ましたが、これは無視して良いらしい。

%% Loading DLL
fname_lib = 'libmecab.dll';
fname_header = 'mecab.h';
[notfound, warnings] = loadlibrary(fname_lib, fname_header);

下のコードで「ライブラリlibmecabdeの関数」がたくさん表示されればDLLを読み込めているということのようです。

%% (optional) view available functions of the DLL
libfunctionsview('libmecab')

あと、下のように2,3行おまじないを書いてから、

%% Calling MeCab morphological analyzer
argv = libpointer('stringPtrPtr', {'MeCab'});
argc = 1;
mecab = calllib('libmecab', 'mecab_new', argc, argv);

日本語のテキストを与えると形態素分析の結果が返ってきます。すばらしい。

%% test MeCab functionality
input = '太郎は次郎が持っている本を花子に渡した。';
result = calllib('libmecab', 'mecab_sparse_tostr', mecab,input);

result =

太郎 名詞,固有名詞,人名,名,*,*,太郎,タロウ,タロー
は 助詞,係助詞,*,*,*,*,は,ハ,ワ
次郎 名詞,固有名詞,人名,名,*,*,次郎,ジロウ,ジロー
が 助詞,格助詞,一般,*,*,*,が,ガ,ガ
持っ 動詞,自立,*,*,五段・タ行,連用タ接続,持つ,モッ,モッ
て 助詞,接続助詞,*,*,*,*,て,テ,テ
いる 動詞,非自立,*,*,一段,基本形,いる,イル,イル
本 名詞,一般,*,*,*,*,本,ホン,ホン
を 助詞,格助詞,一般,*,*,*,を,ヲ,ヲ
花 名詞,一般,*,*,*,*,花,ハナ,ハナ
子 名詞,接尾,一般,*,*,*,子,コ,コ
に 助詞,格助詞,一般,*,*,*,に,ニ,ニ
渡し 動詞,自立,*,*,五段・サ行,連用形,渡す,ワタシ,ワタシ
た 助動詞,*,*,*,特殊・タ,基本形,た,タ,タ
。 記号,句点,*,*,*,*,。,。,。
EOS

最後にDLLを開放します。

%% Release DLL from the memory
clear
unloadlibrary('libmecab')

Matlab (64 bit)でMeCabを使う環境が整ったので、これからいろんなテキストを突っ込んで遊んでみます。