[R言語] Rのお部屋::あーるのおへや[別館]

pythonでは、文字列で与えられた式を評価する場合

eval

を、使うようです．

>>> eval('a')
Traceback (most recent call last):
File "", line 1, in
File "", line 1, in
NameError: name 'a' is not defined
>>> a = 3.5
>>> eval('a')
3.5

式の中で、代入を行う場合はcompileを使うと良いようです．

eval(compile("c = a+b", "(^.^)", "single"))

hogehogeという変数に値を代入する例

>>> eval(compile("c = 1.111", "(^.^)", "single"))
>>> c
1.111
>>> eval(compile("hogehoge = 1.111", "(^.^)", "single"))
>>> hogehoge
1.111
>>> s = "hogehoge = 'xyzzy'"
>>> s
"hogehoge = 'xyzzy'"
>>> eval(compile(s, "(^.^)", "single"))
>>> hogehoge
'xyzzy'

[R] [R言語] [2次計画法] [2次計画問題] [非線形計画法] [数理計画法] [最適化]

Rのお部屋で自分用のメモとして書いておいた記事（２次計画法（２次計画問題）のメモ）を、少しわかりやすくまとめてみることにしました．

自分用のメモのため、見に来てくれた方にはわかりにくいものとなっていましたので．

Rで二次計画法を行うには「quadprog」というライブラリを使います．

quadprogの説明（別ウインドウでPDFを開きます）にノーテーションを合わせます

min(-d^Tb + 1/2 b^TDb)
with constraints A^Tb ≧ b₀

ここで

b ： 求めたい変数のベクトル n×1の行列（要素n個のベクトル）
D ： n×nの対角行列
d ： n×1の行列（要素n個のベクトル）
A ： 制約条件の左辺の係数
b₀ ： 右辺の定数項

1/2 b^TDb

-d^Tb

制約条件は≧としかなっていませんが、等式制約と不等式制約両方扱えます（あたりまえ）．
ここでquadprogの制約条件は行列Aと制約条件の右辺のb₀で設定しますが
等式制約 → 不等式制約
の順でAとb₀を設定しなくてはなりません．これは後で説明します．

関数solve.QPのサンプルの例を実行してみましょう．最初のコメント部分に以下のように書いてあります．

##
## Assume we want to minimize: -(0 5 0) %*% b + 1/2 b^T b
## under the constraints: A^T b >= b0
## with b0 = (-8,2,0)^T
## and
##     (-4 2 0)
## A = (-3 1 -2)
##     ( 0 0 1)
## we can use solve.QP as follows:
##

これはまず求めたい変数の次元が３で、目的関数が

-(0 5 0) %*% b + 1/2 b^T b

制約条件はbを（b1、b2、b3）とすると、 b0 = (-8,2,0)^Tが右辺、Aが左辺ですから

-4×b1 ＋ -3×b2 ＋       ≧ -8
 2×b1 ＋  1×b2 ＋       ≧  2
          -2×b2 ＋ 1×b3 ≧  0

ここで2次の項を構成する対角行列のDがありませんが

1  0  0
0  1  0
0  0  1

Rのソースを記述します．DmatがD行列を表しています．Dは単位行列ですので以下でDmatを設定できます．

Dmat <- matrix(0,3,3)
diag(Dmat) <- 1
Dmat
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

次に目的関数の1次の項のdvecを設定します．ノーテーションではn×1の行列と書きましたが必ず1列なのでRでいうところのただのベクトルで設定します．

dvec <- c(0,5,0)

Amatは上記制約条件のとおりですが、以下の並びになるように設定します．制約条件の不等式の係数の並びとは行と列が逆になっていますので注意が必要です．

     (-4 2 0)
 A = (-3 1 -2)
     ( 0 0 1)

以下でAmatを設定します．

Amat <- matrix(c(-4,-3,0,2,1,0,0,-2,1),3,3)
Amat
     [,1] [,2] [,3]
[1,]   -4    2    0
[2,]   -3    1   -2
[3,]    0    0    1

bvecは制約条件の右辺ですので

bvec <- c(-8,2,0)

実行して

> solve.QP(Dmat,dvec,Amat,bvec=bvec)
$solution
[1] 0.4761905 1.0476190 2.0952381

$value
[1] -2.380952

$unconstrained.solution
[1] 0 5 0

$iterations
[1] 3 0

$Lagrangian
[1] 0.0000000 0.2380952 2.0952381

$iact
[1] 3 2

等式制約がある場合
Rによる二次計画法 例

[R] [R言語] [2次計画法] [2次計画問題] [非線形計画法] [数理計画法] [最適化]

制約条件の順序は、等式制約がある場合前の記事で

等式制約 → 不等式制約

制約条件を以下のように変更してみます（1番目の制約条件を等式制約に変更しました）．

-4×b1 ＋ -3×b2 ＋       ＝ -8
 2×b1 ＋  1×b2 ＋       ≧  2
          -2×b2 ＋ 1×b3 ≧  0

この場合、不等式制約と同じようにRの変数を設定します．

Dmat <- matrix(0,3,3)
diag(Dmat) <- 1
Dmat
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1
dvec <- c(0,5,0)
Amat <- matrix(c(-4,-3,0,2,1,0,0,-2,1),3,3)
Amat
     [,1] [,2] [,3]
[1,]   -4    2    0
[2,]   -3    1   -2
[3,]    0    0    1
bvec <- c(-8,2,0)

実行するときにmeqという引数に等式制約の数を設定します．今回は1個ですので「meq=1」を設定します．

solve.QP(Dmat, dvec, Amat, bvec=bvec, meq=1)
$solution
[1] 1.123596 1.168539 2.337079

$value
[1] -1.797753

$unconstrained.solution
[1] 0 5 0

$iterations
[1] 3 0

$Lagrangian
[1] 0.2808989 0.0000000 2.3370787

$iact
[1] 3 1

これで等式制約を含んだ2次計画法を実行することができました．

-4 × 1.123596 + -3 × 1.168539 -> -8.000001

等式制約条件も満たしています．

Rによる二次計画法 例

[R] [R言語] [2次計画法] [2次計画問題] [非線形計画法] [数理計画法] [最適化]

2次計画法の典型的な例として、3資産を組み合わせてリスクを最小にするようなポートフォリオの各資産のウエイトを求めてみましょう．

３資産を資産A・資産B・資産Cとします．各資産の例えば月次のリターンを

資産A ：  3%、-1.5%、1%、 7%、4%
資産B ： -2%、 1.5%、7%、 8%、3%
資産C ：  4%、 8.5%、3%、-7%、6%

各資産のウエイトですからそれをW_A、W_B、W_Cと置き、また各資産は必ずロングポジションで保有することにすると、制約条件は

W_A＋W_B＋W_C ＝ 1.0
W_A         ≧ 0.0
    W_B     ≧ 0.0
        W_C ≧ 0.0

それではRで各変数を設定していきましょう．
目的関数はW_A、W_B、W_Cで構成されるポートフォリオのリスクですからまずは共分散行列を求めなければなりません．

retA <- c(0.03,-0.015,0.01,0.07,0.04) #資産Aのリターンの系列
retB <- c(-0.02,0.015,0.07,0.08,0.03) #資産Bのリターンの系列
retC <- c(0.04,0.085,0.03,-0.07,0.06) #資産Cのリターンの系列
ret <- cbind(retA,retB,retC)
C <- cov(ret)
C
           retA       retB       retC
retA  0.0010200  0.0004875 -0.0015475
retB  0.0004875  0.0016750 -0.0015750
retC -0.0015475 -0.0015750  0.0035050

Cをそのまま使っても良いのですがDmatに代入しておきます．

Dmat <- C

目的関数に1次の項はありませんので

dvec <- c(0,0,0)

制約条件の左辺から

## W_A＋W_B＋W_C  → (1,1,1)
## W_A          → (1,0,0)
##     W_B      → (0,1,0)
##         W_C  → (0,0,1)

Amat <- cbind(c(1,1,1), c(1,0,0), c(0,1,0), c(0,0,1))
Amat
     [,1] [,2] [,3] [,4]
[1,]    1    1    0    0
[2,]    1    0    1    0
[3,]    1    0    0    1

制約条件の右辺から

bvec <- c(1,0,0,0)

ここで1番目の制約条件だけが等式制約であることを思い出してください．ですのでmeq引数に1をセットしてsolve.QPを実行します．

res <- solve.QP(Dmat, dvec, Amat, bvec, meq=1, factorized=FALSE)

結果を見てみます．

res
$solution
[1] 0.4624505 0.2148457 0.3227038

$value
[1] 3.852634e-05

$unconstrained.solution
[1] 0 0 0

$iterations
[1] 2 0

$Lagrangian
[1] 7.705267e-05 0.000000e+00 0.000000e+00 0.000000e+00

$iact
[1] 1

$solutionが求めたい変数の結果ですので、これがウエイトになっています．合計すると１になっていて制約条件をちゃんと満たしていますね．

sum(res$solution)
[1] 1

これで、これらの資産のリスクが最小になるウエイトを求めることができました．

ここ数年「ビッグデータ」という言葉をよく耳にするようになってきました．

ビッグデータとはなんなのでしょうか．

例えばYahooのページが閲覧されたという閲覧ログ．毎秒数千とか数万というオーダーのPVなのだと思います．もっと多いのかも．

例えば、ソーシャルゲーム．何万という個人がユーザーとして存在し、そのうち誰がどういったボタンを押しただとかそういったページ遷移のログ、あるいは行動履歴．

鉄道会社の各乗客の行動履歴．最近は電子マネーでサービスを利用することが多いので、ある人がいつどの駅で乗ってどの駅で降りたということのログ．

電話会社の通話履歴．ある人がどこから、誰へ（どこへ）、いつ、どれくらいの長さの通話をしたかというログ．

POSデータ．全国のコンビニで毎日たくさんの個人が買い物をします．これも電子マネーを利用する人が多いので誰がいつどのような商品を購入したかというログ．

例えばAmazon．どのユーザーがどんな商品を見たか．そしてどういう商品を見た後に、どういう商品を購入したかというログ．

1秒に数千とか数万というレベルで発生するログは、これまでは残すにしても膨大なストレージが必要でした．またそれを集計しようにも高性能なコンピュータが必要でした．ところがストレージは高性能・低価格化が進みそういったログを残すことが可能になってきました．また膨大なログでもそれを集計するためのハードウエアや仕組みが整備されてきて、集計することもそんなに大変ではなくなってきました．

ビッグデータとはこういった、これまでは残すにしてもコストがかかりすぎるし、集計しようにもログの発生スピードのほうが速く、集計・分析ができないとして、捨てられてきたデータと考えてよさそうです．

こういった膨大なログだとか行動履歴を集計し、分析することでそこから何かビジネスに活かせる知見を得たり、サービスの向上に貢献させるということを具現化させていた例としてあまりにも有名なのがAmazonでしょう．

Amazonのレコメンドは10年くらい前からあったように思います．今でこそ、そのレコメンドには心くすぐられるものがありますが、初めのころはひどいレコメンドもありました．

Amazon風「おすすめURL」にも書きましたので引用します．

最近のＡｍａｚｏｎのおすすめはかなり精度がよいと感じますが、一昔前は、「はぁ？」と感じることも少なくありませんでした．実際にあったことなのですが、ＰＣ用の液晶モニターをＡｍａｚｏｎで注文し購入したことがあったのですが、その後のおすすめに、１８禁のゲームソフトがおすすめされたことがありました．

筆者と同じ液晶モニターを購入した誰かが合わせてそういうのも購入したのだろうと思いますが、さすがにこの「こんなアイテムも購入しています」にはドン引きました．最近ではこういうことは起こらず、同じようなカテゴリーの中からおすすめが表示されているように思います．

現在のＡｍａｚｏｎのレコメンドシステムはかなり工夫が凝らされているようです．

今ではこんなことはないと思います．仮にあっても「ご愛嬌」と思えるくらい、Amazonのレコメンドはついクリックしたくなります．

つまり、世の中はこう考えたのです．

「こういった購入履歴や閲覧履歴を残して、活かすことができるんだ」

Amazonが失敗していたら、ビッグデータというイノベーションは起こらなかったかもしれません．いえそれでも別のところが成功していたかもしれませんが．

もちろん、これまで捨てられていたデータを残したからといって、それだけでなにかビジネスに直結するわけではありません．

pythonで10進数から36進数に基数変換しなければならなかったのでGoogle先生に

「python base36」

で、教えてもらいました．

Python base 36 encoding

いちばん最後のでよいのかな

import string, math

int2base = lambda a, b: ''.join(
    [(string.digits + string.lowercase + string.uppercase)[(a/b**i)%b]
     for i in xrange(int(math.log(a, b)), -1, -1)]
    )

num = 1412823931503067241
test = int2base(num, 36)
test2 = int(test, 36)
print test2 == num

int2baseという関数がそれです．
int2base(x, base)というインターフェースになっています．

確認

>>> int2base(36, 36)
'10'
>>> int2base(1, 36)
'1'
>>> int2base(2, 36)
'2'
>>> int2base(36*36, 36)
'100'

[DaReCo] [連想計算]

DaReCoというウエブアプリケーションを作ってみました．

作った動機といいますか、作るまでの経緯は以下のとおりです．

趣味（？）で連想計算エンジンをpythonで実装していました．メモリベースのものやNoSQL（MongoDBやCassandra）ベースのもの、ファイルシステムを永続化機構として利用したものなど、です．

メモリベースの連想計算エンジンに対しそのパフォーマンスと計算の結果が（自分の）感覚に合うかということを検証するために、いろいろなポータルサイトから記事を集めそれを連想計算エンジンにぶちこんでいました．

ここで連想計算エンジンに入力したのが、「A：記事のURL」と「B：記事のタイトルを形態素解析しその結果から連想されるワード」です．連想計算エンジン自体の計算の確認では、ワードのリストを与え（B側）、そのワードのリストから計算される記事を表示（A側）するということでした．

結果としてある一定の条件下では、計算スピードは相当速かったです．メモリベースなのでスピードが出なければ困るのですが．また計算結果もそんなに悪くありませんでした（自分の感覚とたいして違わないという程度ですが・・・）．

確認は「ビッグデータ」とか「python」とか自分の興味のある言葉でやっていたのですが、出力される結果（記事のリスト）がけっこう自分の興味に刺さるような結果で、「この記事あとで読まなきゃ」と思うことが多かったのです．

ここまではターミナル上での確認ですから、自分の興味のある記事で、今日更新された記事にどういうものがあるかを見るためにはわざわざ自宅のターミナルでパフォーマンス確認用のスクリプトをたたかなければなりません．

それだったら、どこでもチェックできるようにウエブアプリにしてしまおうと、考えたわけです．

Googleで同じようなことができるかと思ったのですが、興味のあるワードで検索しても最新の記事を表示してはくれません．例えば「python」と入力してもpythonの公式サイトとかドキュメントページしか上位には表示してくれませんよね．

Googleアラートなら同じような最新のニュースを教えてくれるかもしれませんが、連想計算ではないのでpythonというワードに引っかかるページは教えてくれますが、pythonという空間に近い空間の記事は表示してくれません．

そんな感じでGoogleだとできないなぁと感じたのも、自分でアプリにしてみようと思った最大の要因でしょうか．

単純にpythonと入力して連想計算させて記事を出力させてもpythonというワードがタイトルに含まれるものだけを出力するということではありません．pythonが属する意味的？な空間に属する記事も出力します．pythonですからLL（Lightweight Language、軽量プログラミング言語）周りの記事も出力してくれるのです．たとえばRubyとかPerlとかjavascriptとかあるいはそういった言語に関する記事を．

連想計算ですのでpythonというワードが含まれる記事が上位に来ますが、pythonというワードがタイトルに含まれない記事でも、pythonという意味空間に属するあるいは近い記事は出力してくれるのです．

そういったLL周りの情報・動向も知っておきたいでしすしね．

この感覚が個人的にとても気に入っています．

とにかく自分がほしいと思った機能をもつアプリケーションを作ってみたら、個人的にとても楽しく活用できるページができてしまった、そんな感じです．

[DaReCo] [連想計算]

[DaReCo] DaReCoというウエブアプリケーションを作ってみましたの続きです．

先の記事には、「個人的にとても気に入っていて、楽しく活用できるアプリができた」と、書きましたが、わたしも人間ですので、こういったアプリが他人からどう見られるのかなとか、連想計算に同じように興味を持つ人と情報交換してみたいなとも、思うわけです．

さて、DaReCoのページを見ていただくとわかりますが、入力したフレーズ（入力は単語でなくても良い．文章でも良いです）から関連ワードを計算し表示しています．これも連想計算です．この関連ワードには表示はされていませんが、重み（ウエイト）を内部では持っています．このワードとそれに紐付くウエイトの組み合わせによって連想計算を行っていますので、入力したワードが含まれていない場合でも近しい空間の記事を表示してくれますし（正確には近しい空間の記事を表示しているように感じるだけ・・・）、入力したワードをタイトルに含むものが記事の上位に来ます．

この関連ワードはWikipedia様の情報を利用していますので、日本語に関してはほぼ問題がない程度といえます．表記ゆれは勘弁してほしいと思います．特定のニッチなサブカルは捉え切れていないのかもしれません．主にWikipediaが・・・

仮に連想ワードが表示されていたとしても、それに関する記事が必ず表示されるわけではありません．だって自分の興味のあるポータルサイトの記事しか集めていませんから．それに記事はRSSフィードで配信されているものだけですから、RSSで配信していないサイトの記事は入っていません．

つまり、わたしの独断と偏見でポータルサイトを選び、またRSSで配信されているもののみという制限があります．あしからず・・・

RSSで取得した記事の保持期間は現在10日です．前の記事に「趣味でやっている」と書きましたが、連想計算は自宅のマシンで運用しています．現在連想計算を行うモジュールは完全メモリベースですので、計算マシンの性能（主にメモリ容量）に依存します．スワップさせてしまうと計算の性能が出ませんので、スワップしない程度の記事数にしておかなければなりません．現在は10日分を保持し、publishedから10日経過した記事は削除しています．

メモリは32Gありますので、まだまだ余裕はありますが、記事を取得するサイトを増やすこともあるかもしれません．また個人的には最新の記事のみをチェックしたいので、10日で十分なのです．毎日どんな記事がリリースされているかチェックしていますので、本当は5日でも良いくらいです．メモリに余裕があるので10日にしています．

現在連想計算は記事のタイトルのみで行っています．本文で行うこともできますが、RSSをソースとしていますので現時点ではすべての記事にその対応を入れることは不可能です．記事には記事の内容にふさわしいタイトルがつけられていることを前提にしています．

前の記事でも述べましたが、ほんと、個人的にはとても満足しています．