Search on the blog

2016年8月28日日曜日

信頼区間付きの線形回帰予測

 pythonで線形回帰の信頼区間付きpredictionをしたかったけど、sklearnには該当する機能が無かった。Gradient Boosting Regressionだと出来るのだけど[2]、単純なLeast Squareでやりたかったのでどうやるのか調べていたら、stasmodelsを使えばいい[1]らしい。

サンプルソース

import numpy as np
import statsmodels.api as sm
from statsmodels.sandbox.regression.predstd import wls_prediction_std
import matplotlib.pyplot as plt

np.random.seed(1234)
nsample = 30
x1 = np.linspace(0, 20, nsample)
x = np.column_stack((x1, x1**0.5))
x = sm.add_constant(x)
beta = [1, 0.5, 0.5]
y = np.dot(x, beta) + np.random.normal(size=nsample)
res = sm.OLS(y, x).fit()
prstd, iv_l, iv_u = wls_prediction_std(res)

plt.plot(x1, y, 'o', label="Data")
plt.plot(x1, res.fittedvalues, 'r--')
plt.plot(x1, iv_u, 'r--', label="Ordinary Least Squares")
plt.plot(x1, iv_l, 'r--')
plt.legend(loc="best");
plt.show()

wls_prediction_stdは、学習データに対する予測値の有意水準5%信頼区間を返してくれる。有意水準は引数alphaで自由に設定できる。また、学習データ以外のデータに対する予測を行いたい場合は、引数exogにデータを渡せばよい[3]。

結果

参考
[1] https://www.reddit.com/r/MachineLearning/comments/3raivl/code_to_calculate_confidence_interval_for_linear/
[2] http://scikit-learn.org/stable/auto_examples/ensemble/plot_gradient_boosting_quantile.html
[3] https://github.com/statsmodels/statsmodels/blob/master/statsmodels/sandbox/regression/predstd.py

2016年8月23日火曜日

INETドメインとUNIXドメイン

ソケットには、INETドメインソケットとUNIXドメインソケットが存在する。

INETドメイン
異なるマシンで動作しているプロセス間の通信を行うためのソケット。

UNIXドメイン
同じマシン内で動作しているプロセス間の通信を行うためのソケット。
INETドメインソケットでもループバックアドレスを利用することで同一マシン内の通信はできるが、UNIXドメインソケットの方が高速である。

INETドメインソケットのサンプル

inet_server.py

bindするときにIP、PORTを指定する。

import socket

IP = '127.0.0.1'
PORT = 5005
BUFFER_SIZE = 1024

s = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
s.bind((IP, PORT))
s.listen(1)
conn, addr = s.accept()
while 1:
 data = conn.recv(BUFFER_SIZE)
 if not data:
  break
 msg = data.decode('utf-8')
 msg = msg[::-1]
 conn.send(msg.encode('utf-8'))
conn.close() 

inet_client.py

import socket

IP = '127.0.0.1'
PORT = 5005
BUFFER_SIZE = 1024

s = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
s.connect((IP, PORT))
msg = "ABCDEFG"
s.send(msg.encode('utf-8'))
data = s.recv(BUFFER_SIZE)
print (data.decode('utf-8'))
s.close()

プロセスIDからINETドメインソケットを調べる。
$ ps -ef | grep python
  501  9165  2011   0 11:51PM ttys000    0:00.04 python inet_server.py
$ lsof -n -i -P | grep 9165
python3.5  9165 kenjih    3u  IPv4 0xf27949601620cf2b      0t0  TCP 127.0.0.1:5005 (LISTEN)

UNIXドメインソケットのサンプル

unix_server.py

bindするときにファイル名を指定する。
UNIXドメインソケットでは、アドレス・名前空間としてファイルシステムを使用する。

import socket
import os

PATH = '/tmp/sample.sock'
BUFFER_SIZE = 1024

s = socket.socket(socket.AF_UNIX, socket.SOCK_STREAM)
os.unlink(PATH)
s.bind(PATH)
s.listen(1)
conn, addr = s.accept()
while 1:
 data = conn.recv(BUFFER_SIZE)
 if not data:
  break
 msg = data.decode('utf-8')
 msg = msg[::-1]
 conn.send(msg.encode('utf-8'))
conn.close() 

unix_client.py

import socket

PATH = '/tmp/sample.sock'
BUFFER_SIZE = 1024

s = socket.socket(socket.AF_UNIX, socket.SOCK_STREAM)
s.connect(PATH)
msg = "ABCDEFG"
s.send(msg.encode('utf-8'))
data = s.recv(BUFFER_SIZE)
print (data.decode('utf-8'))
s.close()

プロセスIDからUNIXドメインソケットを調べる。
$ ps -ef | grep python     
  501 11543  2011   0 12:02AM ttys000    0:00.04 python unix_server.py
$ lsof -n -P -U | grep 11543
python3.5 11543 kenjih    3u  unix 0xf279496017547fcb      0t0      /tmp/sample.sock

2016年8月18日木曜日

httpbinでテストする

Pythonのurllib3のRetry、Timeoutの動作確認をするためにテスト用のサーバを書こうかどうか迷っていたところ、httpbinという便利なサイトを見つけた。
これは便利。

指定したステータスコードを返す
https://httpbin.org/status/:code
をたたくと、指定した:codeをステータスコードとして返してくれます。

urllib3でステータスコードが400、500番台の場合はバックオフしてリトライするコードをhttpbinで動作確認。
from urllib3.util import Retry
from urllib3 import PoolManager

http = PoolManager(retries=Retry(total=4, backoff_factor=1.0, status_forcelist=range(400, 600)))

response = http.request('GET', 'https://httpbin.org/status/500')
print ("status code=%s." % response.status)

指定した秒だけ待ってレスポンスを返す
https://httpbin.org/delay/:n
をたたくと、min(n, 10)秒間待ってレスポンスを返してくれます。

urllib3のTimeoutの動作確認用ソースです。5秒たつとtimeoutエラーが発生します。
from urllib3.util import Retry, Timeout
from urllib3 import PoolManager

http = PoolManager(timeout=Timeout(5))

response = http.request('GET', 'https://httpbin.org/delay/8', retries=False)
print ("status code=%s." % response.status)

timeoutしても頑張るバージョン。
from urllib3.util import Retry, Timeout
from urllib3 import PoolManager

timeout = Timeout(5)
retries = Retry(total=4, backoff_factor=1.0, status_forcelist=range(400, 600))

http = PoolManager(timeout=timeout, retries=retries) 
response = http.request('GET', 'https://httpbin.org/delay/8')
print ("status code=%s." % response.status)
この場合は、4回リトライしtimeoutエラーになります。

2016年8月15日月曜日

多項ロジットの理論と実装

問題設定
訓練データ(xn, yn), x∈ R D, y∈ {1,2,...C}, n=1,2,...,Nを使って、マルチクラス識別問題を解くモデルを考えます。

多項ロジットのモデル
多項ロジットでは以下のモデルを使います。


分子は訓練データxの特徴量の線形結合をロジスティック関数に入れたものです。
分子は非負の値を取ります。

上の式はxがクラスyに属する確率を表します。
よって、xがクラス1に属する確率、xがクラス2に属する確率、...をすべて足すと1.0になっていて欲しいです。この条件を満たすように分母の値で割って正規化を行います。

パラメータの学習
多項ロジットのパラメータは、


です。ここでCはクラスの個数、Dはデータの次元数です。
つまりクラスの数だけ重みベクトルが存在します。

これらのパラメータを学習するために、以下の対数尤度関数を最大化します。

対数尤度関数の偏微分を計算します。
まず第一項の偏微分は以下のようになります。

[]は中の条件が真のときは1、偽のときは0になるような関数です。

第二項の偏微分は以下のようになります。


よって全体の偏微分は、
となります。

偏微分が求まったので、θを適当な値に初期化して、勾配法で対数尤度関数を最適化することができます。

実装
Pythonで実装しました。

import numpy as np

class MultinomialLogit:

 def __init__(self):
  self.C = 0
  self.N = 0
  self.D = 0
  self.theta = None

 def prob(self, X, y):
  return np.dot(self.theta[y], X) / sum(np.dot(self.theta[c], X) for c in range(self.C))

 def logLiklihood(self, X, y):
  return sum(np.log(self.prob(X[n], y[n])) for n in range(self.N))
 
 def fit(self, X, y, iter=200):
  self.C = len(set(y))
  self.N = len(X)
  self.D = len(X[0])+1   # bias term is added
  self.theta = np.random.rand(self.C, self.D);

  bias = np.ones((len(X), 1))
  X = np.hstack((X, bias))

  def delta(a, b):
   return 1.0 if a == b else 0.0
  
  for i in range(iter):
   for c in range(self.C):
    gradient = sum((delta(c, y[n]) - self.prob(X[n], c)) * X[n] for n in range(self.N)) 
    self.theta[c] += gradient
   print ("i=%d %f" % (i, self.logLiklihood(X, y)))
    
 def predict(self, X):
  bias = np.ones((len(X), 1))
  X = np.hstack((X, bias))
  return [np.argmax([self.prob(x, c) for c in range(self.C)]) for x in X]

 def score(self, X, y):
  Ypred = self.predict(X)
  comp = Ypred == y
  return sum(comp) / len(comp)

Irisで実験
Irisのデータを使って実験してみました。
テスト識別率は80%後半程度です。sklearnのLogisticRegressionだと100%が出るので全然ダメな感じですが、
  • iterationごとに対数尤度関数が上昇すること
  • 対数尤度関数が上がるとテスト識別率も上がること
を確認できたのでそれっぽい実装は出来てそうです。学習率を設定したり、 ヘッセ行列を計算してニュートン法系の最適化を行ったり、確率的勾配法を使ったりなど改善の余地はたくさんありそうです。

if __name__ == '__main__':
 from sklearn import datasets
 from sklearn import cross_validation
 
 iris = datasets.load_iris()
 X = iris.data
 y = iris.target

 Xtrain, Xtest, Ytrain, Ytest = cross_validation.train_test_split(X, y, test_size=0.3)
 logit = MultinomialLogit()
 logit.fit(Xtrain, Ytrain, iter=1000)
 print (logit.score(Xtest, Ytest))


2016年8月14日日曜日

SRM 507 Div1 500 CubePacking

問題
辺の長さが1の立方体がNs個、辺の長さがLの立方体がNb個ある。
これらの立方体を出来るだけ体積が小さい直方体に格納したい。
直方体の体積の最小値を求めよ。

解法
解の下限はNb*L*L*L+Nsである。
解の上限はNb*L*L*L+Ns+L*L-1である。(L*Lの底面を持つ直方体に立方体を詰め込む)
ということで、 解の候補をすべて列挙しても高々L*L程度のループをまわせばよい。

あとは、それぞれの解の候補に対して、その体積になるような直方体を全列挙すればよい。これは時間がかかりそうだが約数の個数は思ったほど多くないため十分に間に合う。

以下に整数の約数の個数がどの程度かを示す。
この結果はN周辺の数をサンプリングし、それらの約数の個数の平均値を取ったものである。
N 約数の個数
106 14.9
107 17.3
108 19.6
109 21.9

ざっくりと高々log(N)程度と覚えておくと良さそう(あくまでも平均値であることに注意)。

ソースコード
using namespace std;

#define ALL(x) (x).begin(), (x).end()
#define EACH(itr,c) for(__typeof((c).begin()) itr=(c).begin(); itr!=(c).end(); itr++)  
#define FOR(i,b,e) for (int i=(int)(b); i<(int)(e); i++)
#define MP(x,y) make_pair(x,y)
#define REP(i,n) for(int i=0; i<(int)(n); i++)

vector<long long> divisor(long long n) {
  vector<long long> ret;
  for (long long i = 1; i * i <= n; i++) {
    if (n % i == 0) {
      ret.push_back(i);
      if (i != n / i)
        ret.push_back(n / i);
    }
  }
  return ret;
}

class CubePacking {
  public:
  int getMinimumVolume(int Ns, int Nb, int L) {
    int s = L*L*L*Nb+Ns;
    for (int x = s;; x++) {
      auto ds = divisor(x);
      REP (i, ds.size()) REP (j, i+1) {
        int a = ds[i];
        int b = ds[j];
        int c = x / a / b;
        if (a * b * c != x)
          continue;
        if ((a/L) * (b/L) * (c/L) >= Nb)
          return x;
      }
    }
    return -1;
  }
};

2016年8月12日金曜日

TornadoでWebSocketを使ってみた

WebSocketとは?
ウェブサーバとウェブブラウザ間で双方向通信を行うための規格。
WebSocketを使うことで、サーバ側からクライアント側にプッシュ配信を行うことができる。通信ごとにコネクションを張り替えるのではなく、すべての通信を一つのコネクションで行うことができる。

Tornadoとは?
PythonのWebフレームワーク/非同期通信のライブラリ。
数万単位のコネクションにスケールすることができ、WebSocketのようにユーザごとに長期のコネクションを構築しておきたい場合に有用。

サンプル
チャットアプリを作ってみた。

サーバ側のコード。
tornado.websocket.WebSocketHandlerクラスを継承して、必要なイベントハンドラメソッドをオーバーライドすればよい。

import tornado.ioloop
import tornado.web
import tornado.websocket

class ChatHandler(tornado.websocket.WebSocketHandler):
 clients = []

 def open(self):
  if self not in ChatHandler.clients:
   ChatHandler.clients.append(self)

 def on_message(self, message):
  for c in ChatHandler.clients:
   c.write_message(message)

 def on_close(self):
  if self in ChatHandler.clients:
   ChatHandler.clients.remove(self)

class MainHandler(tornado.web.RequestHandler):
 def get(self):
  self.render('main.html')

if __name__ == "__main__":
 application = tornado.web.Application([
  (r"/chat/main.html", MainHandler),
  (r"/chat/websocket", ChatHandler),
 ])
 application.listen(8888)
 tornado.ioloop.IOLoop.current().start()

クライアント側のコード。WebSocketクラスを使って接続、メッセージ送信、メッセージ受信時の処理を行う。
<!DOCTYPE html>
<html>
  <head>
    <title>websocket chat</title>
    <script src="http://ajax.googleapis.com/ajax/libs/jquery/1.12.4/jquery.min.js"></script>
    <link rel="stylesheet" href="//maxcdn.bootstrapcdn.com/bootstrap/3.3.2/css/bootstrap.min.css" />
    <script src="//maxcdn.bootstrapcdn.com/bootstrap/3.3.2/js/bootstrap.min.js"></script>
    <script>
      $(function() {
        socket = new WebSocket("ws://kenjih.com:8888/chat/websocket");
        socket.onopen = function() {};
        socket.onmessage = function(e) {
          $("#msgbox").append($('<p>'+e.data+'</p>'));
        };
        $('#send').click(function(){
          socket.send($('#msg').val());
          $('#msg').val('');
        });
        $('#msg').keydown(function(e) {
          if (e.keyCode == 13) {
            socket.send($('#msg').val());
            $('#msg').val('');
          }
        });
      });
    </script>
  </head>
  <body>
    <div class="container">
      <div class="row">
      <h1>Websocket Chat Room</h1>
      </div>
      <div class="row">
        <div class="col-xs-5">
          <div class="input-group">
            <input id="msg" type="text" class="form-control" placeholder="put your message...">
            <span class="input-group-btn">
              <button id="send" class="btn btn-primary" type="button">Send</button>
            </span>
          </div>
        </div>
        <div class="col-xs-7">
          <div class="panel panel-primary">
            <div class="panel-heading">
              <h3 class="panel-title">Messages</h3>
            </div>
            <div id="msgbox" class="panel-body"></div>
          </div>
        </div>
      </div>
    </div>
  </body>
</html>

SRM 502 Div1 500 TheProgrammingContestDivOne

問題概要
プロコンでは問題がn問出題される。
それぞれの問題に対して、
  • point[i] 問題iを解いたときにもらえる得点
  • minus[i] 問題iの単位時間あたりの得点の減衰
  • time[i] 問題iを解くのに必要な時間
が与えられる。問題を解くのに使える時間はTである。
最適な戦略を用いたときに得られる得点の最大値を求めよ。

解法
この問題の面白いところは、
  • 問題を解く順番の最適化
  • 解くべき問題集合の最適化
という2つの最適化を考えなければいけないところ。

もし解く順序が分かっていれば0-1ナップサック問題に帰着できるので、順序を考える。
i=0,1,2,..,n-1の順で解くことを考えてみる。
このとき、
minus[k+1] * time[k] < minus[k] * time[k+1]
を満たさなければ、kとk+1の順序は逆にしたほうがいいことがわかる。

つまり問題を解く最適な順序は、
time[i] / minus[i]
が小さい順ということが分かる。

よって上の順でソートして、0-1ナップサック問題を解けばよい。

ソースコード

#define ALL(x) (x).begin(), (x).end()
#define EACH(itr,c) for(__typeof((c).begin()) itr=(c).begin(); itr!=(c).end(); itr++)  
#define FOR(i,b,e) for (int i=(int)(b); i<(int)(e); i++)
#define MP(x,y) make_pair(x,y)
#define REP(i,n) for(int i=0; i<(int)(n); i++)

long long dp[100000+1];
class TheProgrammingContestDivOne {
public:
  int find(int T, vector<int> point, vector<int> minus, vector<int> time) {
    int n = point.size();

    // optimize order
    vector<pair<double, int> > vs;
    REP (i, n) vs.push_back(make_pair((double)time[i]/minus[i], i));
    sort(vs.begin(), vs.end());
    
    // optimize set
    memset(dp, 0, sizeof(dp));
    for (auto &pr: vs) {
      int i = pr.second;
      for (long long j = T; j >= time[i]; j--) {
        dp[j] = max(dp[j], dp[j-time[i]] + point[i] - minus[i] * j);
      }
    }
    return *max_element(dp, dp+T+1);
  }
};

2016年8月11日木曜日

Apache Thriftを使ってみた

前回のProtocol Buffersに続いて、Thriftを使ってSerialize/Deserialize処理を書いてみた。

インストール
OS Xの場合は、homebrew一発。
$ brew install thrift

サンプルプロジェクトの作成
Javaでサンプルを書くので、gradleプロジェクトを作っておく。

$ gradle init --type java-library

build.gradleに 以下を追記。
apply plugin: 'application'
mainClassName = 'Main'
compile group: 'org.apache.thrift', name: 'libthrift', version: '0.9.3'

IDLファイルの記述
以下のようなthrift IDLを書いて、sample.thriftというファイル名でプロジェクト直下に保存。
namespace java sample

struct Person {
  1: string name,
  2: i32 age,
  3: string sex
}

thriftコンパイラを使って、IDLからソースファイルを自動生成。
thrift --gen java -out src/main/java sample.thrift

ファイルが作成されたことを確認。
$ find . -name "Person.java"
./src/main/java/sample/Person.java

サンプル
サンプル実行ファイルを書く。
import org.apache.thrift.TException;
import org.apache.thrift.protocol.TBinaryProtocol;
import org.apache.thrift.TSerializer;
import org.apache.thrift.TDeserializer;

import sample.Person;

public class Main {
  public static void main(String[] args) throws TException {
    // serialize
    Person p = new Person("taro", 20, "male");
    TSerializer serializer = new TSerializer(new TBinaryProtocol.Factory());
    byte[] bytes = serializer.serialize(p);

    // deserialize
    TDeserializer deserializer = new TDeserializer(new TBinaryProtocol.Factory());
    Person q = new Person();
    deserializer.deserialize(q, bytes);
    System.out.println("Deserialized object: " + q);
  }
}

実行してみる。
$ ./gradlew run
Deserialized object: Person(name:taro, age:20, sex:male)
ちゃんと動いた!

2016年8月9日火曜日

Protocol Buffersを使ってみた

JavaでGoogleのProtocol Buffersを使ってみました。

インストール
以下の公式ページからコンパイラをダウンロード。
https://developers.google.com/protocol-buffers/docs/downloads

READMEに従ってビルド。
$ ./configure
$ make
$ make check
$ make install

IDLからJavaソースを自動生成
作業フォルダに移動して、gradleでプロジェクトを作成する。
gradle init --type java-library

protoファイルを書かないといけないらしいので書く。protoファイルは、いわゆるIDL(インターフェース記述言語)の一種。
$ cat protos/sample.proto 
package sample;

option java_package = "com.kenjih.sample";
option java_outer_classname = "SampleProtos";

message Person {
  required string name = 1;
  required int32 age = 2;
  required string sex = 3;
}

コンパイルして、IDLからjavaファイルを生成する。
$ protoc --java_out=src/main/java protos/sample.proto

なんか出来た。
$ find . -name "*.java"
./src/main/java/com/kenjih/sample/SampleProtos.java

ビルドする。
build.gradleに以下を追記。
apply plugin: 'application'
mainClassName = 'com.kenjih.sample.Main'
dependencies {
    compile group: 'com.google.protobuf', name: 'protobuf-java', version: '3.0.0'
}

./gradlew buildでコンパイル。
BUILD SUCCESSFULとでればOK。

動作確認
動作確認用のコードを書いてみる。
package com.kenjih.sample;

import java.io.FileInputStream;
import java.io.FileOutputStream;
import java.io.IOException;

import com.kenjih.sample.SampleProtos.Person;

public class Main {
  public static void main(String[] args) throws IOException {
    Person.Builder builder = Person.newBuilder();
    Person p = builder.setName("Taro").setAge(20).setSex("Male").build();
    p.writeTo(new FileOutputStream("person.bin"));

    Person q = Person.parseFrom(new FileInputStream("person.bin"));
    System.out.println(q.toString());
  }
}

gradleから実行。
$ ./gradlew run
name: "Taro"
age: 20
sex: "Male"

serialize/deserializeが出来ていることを確認できた。
ついでにバイナリファイルの中身を確認してみる。
$ hexdump -C person.bin
00000000  0a 04 54 61 72 6f 10 14  1a 04 4d 61 6c 65        |..Taro....Male|
0000000e
なるほど。

2016年8月6日土曜日

木の直径を求める

グラフG(V, E)について考える。G(V, E)の頂点間距離の最大値をG(V, E)の直径と呼ぶ。
木の直径はDFSを2回行うことで簡単に計算できる。

AOJに木の直径を求める問題があったので、解いてみた。

using namespace std;

#define REP(i,n) for(int i=0; i<(int)(n); i++)
#define FOR(i,b,e) for (int i=(int)(b); i<(int)(e); i++)
#define ALL(x) (x).begin(), (x).end()

vector<pair<int, int> > g[100000];

pair<int, int> tddfs(int v, int par = -1) {
  pair<int, int> ret = {0, v};
  for (auto &x : g[v]) {
    int w, cost;
    tie(w, cost) = x;
    if (w == par)
      continue;
    auto r = tddfs(w, v);
    ret = max(ret,  {cost+r.first, r.second});
  }
  return ret;
}

int treeDiameter() {
  auto v = tddfs(0);
  auto w = tddfs(v.second);
  return w.first;
}

int main() {
  int n;
  scanf(" %d", &n);
  REP (i, n-1) {
    int x, y, d;
    scanf(" %d %d %d", &x, &y, &d);
    g[x].push_back({y, d});
    g[y].push_back({x, d});
  }
  int ret = treeDiameter();
  printf("%d\n", ret);
  
  return 0;
}