No science, No life.

もっと科学を身近に

相互情報量を計算してみた

「脳の情報表現」という本にある「相互情報量 (mutual information)」を計算してみました。

脳の情報表現―ニューロン・ネットワーク・数理モデル

脳の情報表現―ニューロン・ネットワーク・数理モデル

エントロピーは、ある確率変数 X のもつ不確かさを表す量で、下のように定義されます。

 H(X) = \sum -P(x) \log P(x)

2変数の同時確率の場合は、

 H(X1, X2) = \sum -P(x1, x2) \log P(x1, x2)

です。そして、相互情報量は、

 I(X1, X2) = H(X1) + H(X2) - H(X1, X2)

と定義されます。この相互情報量という値は、一方の確率変数 X1 の値を知ったときに、もうひとつの確率変数 X2 についての不確かさの減少分を指しています。

この本の20ページ目にある下の例について、実際に計算してみました。

X1 と X2 が相関している場合。相互情報量 I(X1,X2) は 1.08 ビットとなりました (片方の X1 の値を知ると、1.08 ビット分だけ、X2 の値の不確かさが減る)。
f:id:neocortex:20190616233827p:plain

X1, X2 の分布は上と同じですが、片方をシャッフルして独立にした場合。相互情報量は低下して、0.10 ビットとなりました (片方の X1 の値を知っても、もう一方の X2 についての情報はほとんど得られない)。
f:id:neocortex:20190616233842p:plain

Matlab コードは下の通りです。

function mi()

xaxis = -3.9:0.2:3.9;
nbin = numel(xaxis);
bin = xaxis(2)-xaxis(1);

% 2d gaussian random variables
X = mvnrnd([0 0],[0.8 0.7; 0.7 0.8],5000);

% shuffling
if 0
    idx = randperm(size(X,1));
    X = [X(:,1) X(idx,2)];
end

% probability distribution function
pdf1 = hist(X(:,1),xaxis); pdf1 = pdf1/sum(pdf1); 
pdf2 = hist(X(:,2),xaxis); pdf2 = pdf2/sum(pdf2);
pdf12 = zeros(nbin,nbin);
for ii=1:nbin
    for jj=1:nbin
        pdf12(ii,jj) = sum(xaxis(ii)-bin/2<X(:,1) & X(:,1)<=xaxis(ii)+bin/2 &...
            xaxis(jj)-bin/2<X(:,2) & X(:,2)<=xaxis(jj)+bin/2);
    end
end
pdf12 = pdf12./sum(sum(pdf12));

% entoropy
H1 = getEntropy(pdf1);
H2 = getEntropy(pdf2);
H12= getEntropy(pdf12);

% mutual information
I12 = H1+H2-H12;

% plot
figure(1); clf
subplot(3,3,[4 8])
scatter(X(:,1),X(:,2),'.')
axis([-4 4 -4 4])
title(sprintf('H(x1,x2) = %.2f, I(x1,x2) = %.2f',H12,I12))
xlabel('x1');
ylabel('x2')

subplot(3,3,[1 2])
bar(xaxis,pdf1,1)
xlim([-4 4])
title(sprintf('H(x1) = %.2f',H1))
ylabel('P(x1)')

subplot(3,3,[6 9])
barh(xaxis,pdf2,1)
ylim([-4 4])
title(sprintf('H(x2) = %.2f',H2))
xlabel('P(x2)')

end

function H = getEntropy(pdf)

% linearlize
pdf = pdf(:);

% normalize (just in case)
pdf = pdf./sum(pdf);

% remove 0
pdf(pdf==0) = [];

% entropy 
H = sum(-pdf.*log2(pdf));

end

MatlabでTwitterのAPIを使う (3)

Matlabツイッターから情報を集めるには、twittyがとても便利です。ただ、つぶやきを投稿するのには、自分の環境ではエラーが出てうまく使えませんでした。

さがしてみると、Update twitter statusというそのものズバリなコードがあり、これで上手くツイートすることができました。
jp.mathworks.com
ファイルをダウンロードしたら、まず、twit.mを 1 行だけ書き換えます。下の行を探して、.../1/を.../1.1/としましょう。

% url = URL('https://api.twitter.com/1/statuses/update.json');
url = URL('https://api.twitter.com/1.1/statuses/update.json');

それから、下のように入力するとGUIが出てくるので、ツイッターのアクセストークンを入力・保存します。

>> twitpref

それから、下のコードでツイートできます。すばらしい。

>> twit('テスト');

データをクラスタリングするk-means法を書いてみた

データをいくつかのグループにクラスタリングする手法としてk-means法があります。MatlabのStatistics and Machine Learning Toolboxには、kmeansというそのものズバリの関数があるのですが、アルゴリズムを体感するために自分で書いてみました。

Wikipediaによると、大まかなアルゴリズムは下の通り。

  1. 各データ点x(i)に対してランダムにクラスタを割りあてる。(クラスタ数は最初に自分で決める必要がある)
  2. 割りあてたデータ点をもとに、各クラスタの中心V(j)を計算。
  3. 各x(i)と各V(j)との距離を求め、x(i)を最も近い中心のクラスタに割りあてなおす。
  4. 上記の処理で、すべてのx(i)のクラスタの割り当てが変化しなかった場合、収束したと判断して終了。そうでない場合は、新しく割り振られたクラスタからV(j)を再計算して上記の処理を繰り返す。

これをそのまま試してみました。データとしては、50点からなるクラスタを適当に4つ撒きました。

f:id:neocortex:20170612005002p:plain

f:id:neocortex:20170612005014p:plain

f:id:neocortex:20170612005027p:plain

f:id:neocortex:20170612005038p:plain

f:id:neocortex:20170612005050p:plain

f:id:neocortex:20170612005101p:plain

f:id:neocortex:20170612005118p:plain

このデータだと7ステップで収束して処理が終了しました。それっぽくクラスタリングできていますね。

このオリジナルの方法はとても単純なのですが、このままだと、一番最初のランダムなクラスタの割りあてに依存して、明らかに変なクラスタリングをしてしまう場合もあります。この点を改良したk-means++法というのもあり、Matlabのkmeans関数で実装されているのはこちらのようです。

ほかに、複数のクラスタへの帰属を許した fuzzy c-means法というのもあったりします。

今回のコードは下の通り。

clear
rng(13)

% data
x = randn(50,2);
for ii=1:3
    x = [x; randn(50,2)+repmat(rand(1,2)*5,50,1)];
end

% number of clusters 
k = 4;

% assign initial clusters randomly
clu = randi(k,[size(x,1),1]);

h = figure(1); clf
cycle=0;
while 1
    cycle = cycle+1;
    
    % centroids of clusters
    V = nan(k,2);
    for ii=1:k
        V(ii,:) = mean(x(clu==ii,:));
    end
    
    % plot 
    for ii=1:k
        rgb = hsv2rgb(1/k*ii,1,1);
        % individual data points
        plot(x(clu==ii,1),x(clu==ii,2),'ko','markerfacecolor',rgb); hold on
        % centroid
        plot(V(ii,1),V(ii,2),'*','markersize',15,'color',rgb);
    end
    title(cycle)
    drawnow
    hold off
    % save image
    saveas(h,num2str(cycle),'png')
    
    % distance between data points and centroids
    d = nan(size(x,1),k);
    for ii=1:k
        d(:,ii) = sqrt( (x(:,1)-V(ii,1)).^2+(x(:,2)-V(ii,2)).^2 );
    end
    
    % update cluster label to the nearest centroid
    [~,cluNew] = min(d,[],2);
    
    % break if the new cluster label is identical to the previous one
    if any(clu~=cluNew)
        clu = cluNew;
    else 
        break;
    end
end

自己組織化マップ (self-organizing map) を書いてみた

情報を自動で分類するアルゴリズムとして、自己組織化マップ (self-organizing map) というものが知られています。自己組織化マップは教師無し学習の一種で、おおざっぱな学習則はつぎの通りです。

1. ランダムな重みベクトルからなるマップを用意する。
2. 入力ベクトルを1つ用意する。
3. 入力にもっとも近いノードをみつける。
4. そのノードと、近傍のノードとを、入力ベクトルに近づける。
5. ステップ2-4をたくさん繰り返す。

ちゃんとした学習則は、Wikipediaなどを参考に。提唱者Kohonen自身による下の著作もあります。

自己組織化マップ

自己組織化マップ

これをMatlabで書いてみました。色 (赤・緑・青の3要素からなるベクトル) を自動で分類します。最初のランダムなマップはたとえば下のような状態。

f:id:neocortex:20170212233639p:plain:w300

ここから、ステップ2-4を繰り返していきます。100回目が終わるとこんな感じ↓ まだあまり変化が分かりません。

f:id:neocortex:20170212233815p:plain:w300

200回目↓

f:id:neocortex:20170212234355p:plain:w300

500回目↓ マップがすこし滑らかになってきました。

f:id:neocortex:20170212234045p:plain:w300

1000回目↓

f:id:neocortex:20170212234129p:plain:w300

2000回目↓

f:id:neocortex:20170212235255p:plain:w300

5000回目↓ 似た色が自動的に近くのノードに配置されています。

f:id:neocortex:20170212234201p:plain:w300

アニメーションでみるとこんな感じ↓

f:id:neocortex:20170202004023g:plain

シンプルな学習則なのに、情報を自動で分類できるところが面白いですね。「色」は構造のよく分かったデータですが、全体像の見えにくい高次元のデータにも適用することができます (というより自己組織化マップが真価を発揮するのはそういう場面のはず)。

今回のコードは下の通りです (ステップ4の部分はもっとスッキリ書けると思いますが、ひとまずこのままで)。

function SOMtest

figure(1);
clf
N = 20;      % map size
a = 0.1;     % learning coefficient

% random weight vector (step 1)
W = rand(N,N,3);
image(W)
title('0')
axis image off
drawnow
saveas(h,'0','png')

for t=1:5000
    % input vector (step 2)
    V = rand(1,1,3);
    
    % find the nearest node (step 3)
    D = sum((W-repmat(V,N,N,1)).^2,3);
    [~,I] = min(D(:));
    [Irow,Icol] = ind2sub(size(D),I);

    % update the weight vector (step 4)
    for ii=-1:1
        for jj=-1:1
            if 0<Irow+ii && Irow+ii<=N && 0<Icol+jj && Icol+jj<=N
                W(Irow+ii,Icol+jj,:) = ...
                    W(Irow+ii,Icol+jj,:) + a*(V-W(Irow+ii,Icol+jj,:));
            end
        end
    end
    image(W)
    title(t)
    axis image off
    drawnow
    
    % save image
    if mod(t,100)==0
        saveas(h,num2str(t),'png')
    end
end

MatlabでTwitterのAPIを使う (2)

前回の記事のままでは、(英語のツイートは良いのですが) 日本語のツイートはちゃんと表示されません。

たとえば下のように@MATLAB_jpのツイートを取得すると、\u(4ケタの16進数)の羅列からなるユニコードが得られます。

tw = twitty();
S = tw.userTimeline('screen_name', 'MATLAB_jp');
disp(S{1}{1}.text)

#MATLAB \u4fbf\u5229\u95a2\u6570\u7d39\u4ecb \u300cismember\u300d...

そこで、下の関数unicode2txt.mで読めるテキストに変換することにします。関数hex2decでユニコードの16進数を10進数に変換し、さらに関数charで文字に変換しています。

% convert unicode-including string to readable text
function str = unicode2txt(str)

idx = strfind(str,'\u');
for ii=1:length(idx)
    idxTmp = idx(end+1-ii);
    str(idxTmp) = char(hex2dec(str(idxTmp+2:idxTmp+5)));
    str(idxTmp+1:idxTmp+5) = [];
end

これで下のような日本語テキストを得ることができます。

tw = twitty();
S = tw.userTimeline('screen_name', 'MATLAB_jp');
disp(unicode2txt((S{1}{1}.text)))

#MATLAB 便利関数紹介 「ismember」 はデータの重複した要素を探す関数ですが、...

MatlabでTwitterのAPIを使う

MatlabTwitter APIを使ってデータを取得する方法をまとめました。

Twitterのアカウントを持っていない場合は、まずアカウントを作ります。それから、Twitter Platformにアクセスして、APIを使うために必要なキー4種類を取得します。このサイトのお世話になりました。

次に、File ExchangeからtwittyJSON Parserをダウンロードします。このtwittyがスグレモノです。parse_json.mはtwittyのファイル群と同じフォルダに置いておきましょう。

で、Matlabを起動して、twitty.mの説明通りに設定します。まずは、

>> tw = twitty()

と入力します。何か警告が出るかもしれませんが、続けて

>> tw.saveCredentials()

と入力。すると下のようなGUIが登場するので、先ほど取得したキーをコピペします。このキーの入力は1回限りです (Matlabを再起動しても再び入力する必要はない)。

f:id:neocortex:20161228235942p:plain

これで準備完了です!試しに、@MATLABのツイートを取得してみましょう。

>> S = tw.userTimeline('screen_name', 'matlab')

結果は下のようになって、S{1}{インデックス}.textにツイートの内容が格納されます。

>> S

S = 

    {1x20 cell}

>> S{1}{1}

ans = 

                   created_at: 'Wed Dec 28 13:24:16 +0000 2016'
                           id: 8.1410e+17
                       id_str: '814099651293020160'
                         text: 'Get started with MATLAB and #Ardui...'
                    truncated: 0
                     entities: [1x1 struct]
            extended_entities: [1x1 struct]
                       source: '\u003ca href="http://www.hootsuite...'
        in_reply_to_status_id: []
    in_reply_to_status_id_str: []
          in_reply_to_user_id: []
      in_reply_to_user_id_str: []
      in_reply_to_screen_name: []
                         user: [1x1 struct]
                          geo: []
                  coordinates: []
                        place: []
                 contributors: []
              is_quote_status: 0
                retweet_count: 11
               favorite_count: 23
                    favorited: 0
                    retweeted: 0
           possibly_sensitive: 0
                         lang: 'en'

おもなメソッドは、twittyのページに書いてありますし、一覧は下のようにタイプすると得られます。簡単ですね。

>> twitty.API
%% METHODS INTERFACING TWITTER API:            
%% TIMELINES
 homeTimeline()              Returns a collection of the most recent statuses posted by the authenticating user 
 userTimeline()              Returns a collection of the most recent Tweets posted by the user indicated by the
 mentionsTimeline()          Returns the 20 most recent mentions (status containing @username) 
 retweetsOfMe()              Returns the most recent tweets of the authenticated user that have been retweeted
%% STREAMING
 sampleStatuses()            Returns a small random sample of all public statuses via streaming API.
 filterStatuses()            Returns public statuses that match one or more filter predicates.
%% TWEETS
 updateStatus()              Update the authenticating user's status (twit). 
 retweets()                  Returns up to 100 of the first retweets of a given tweet.
 showStatus()                Returns a single status, specified by the id parameter below.
 destroyStatus()             Destroys the status specified by the required ID parameter.
 retweetStatus()             Retweets a tweet specified by the required ID parameter.
 retweeters()                Returns a collection of up to 100 ids of users who retweeted the specified status.
%% SEARCH
 search()                    Search twitter.
%% FRIENDS & FOLLOWERS
 friendsIds()                Returns an array of numeric IDs for every user the specified user is following.
 followersIds()              Returns an array of numeric IDs for every user following the specified user.
 friendshipsLookup()         Returns the relationship of the authenticating user to a list of up to 100
 friendshipsCreate()         Allows the authenticating users to follow the specified user.
 friendshipsDestroy()        Allows the authenticating users to unfollow the specified user.
 friendshipsShow()           Returns detailed information about the relationship between two arbitrary users.
%% USERS
 usersLookup()               Return up to 100 users worth of extended information, specified by either ID, 
 usersShow()                 Returns extended information about a given user.
 usersSearch()               Runs a search for users similar to "Find People" button on Twitter.com.
%% ACCOUNTS
 accountSettings()           Returns settings (including current trend, geo and sleep time information) 
 accountVerifyCredentials()  Test if supplied user credentials are valid.
%% TRENDS
 trendsPlace()               Returns the top 10 trending topics for a specific place (woeid), 
 trendsAvailable()           Returns the locations that Twitter has trending topic information for.
 trendsClosest()             Returns the locations that Twitter has trending topic information for, closest to
%% PLACES & GEO
 geoInfo()                   Returns all the information about a known place.
 geoReverseCode()            Given geographical coordinates, searches for up to 20 places that can be used 
 geoSearch()                 Search for places that can be attached to a statuses/update.
 geoSimilarPlaces()          Locates places near the given coordinates which are similar in name.
%% HELP
 helpConfiguration()         Returns the current configuration used by Twitter including twitter.com slugs 
 helpLanguages()             Returns the list of languages supported by Twitter along with their ISO 639-1 code.
%% THE MAIN API CALLER
 callTwitterAPI()            Call to the twitter API.

カオス (2)

Pythonの勉強のために、前回にMatlabで書いたロジスティック写像スクリプトを、Python 3.5で書き直してみます。Matlabの1から始まるインデックスに慣れているので、Pythonの0から始まるインデックスは混乱します・・・。

import matplotlib.pyplot as plt

# 漸化式の計算
a = 4
x = [0.9]
for i in range(99):
    x.append(a * x[i] * (1-x[i]))

# プロット
plt.plot(range(1,101,1),x)
plt.title('a = '+str(a))
plt.xlabel('n')
plt.ylabel('x[n]')

これで、下のようにa = 4のときの振舞いを描くことができました。

f:id:neocortex:20161204151818p:plain

つぎに、aの値を少しずつ変えたときの分岐図も描いてみます。

import matplotlib.pyplot as plt
import numpy as np

for a in np.arange(0,4,0.001):
    x = [0.9]
    for i in range(99):
        x.append(a * x[i] * (1-x[i]))
    plt.plot([a]*20,x[80:100],'b.')
    
plt.xlabel('a')
plt.ylabel('x[n]')

これで、下のように分岐図をかけました。Pythonはド素人ですが、ちょこちょこ勉強していこうと思います。

f:id:neocortex:20161204160140p:plain

入門 Python 3

入門 Python 3