Movie::くじら前夜

やっとこの記事を書けるくらいまで余裕がでてきました。

えーと、11/30(日)にて行われました、東北大学学友会映画部冬季自主制作上映会「くじら前夜」のオープニング映像を作成しましたので、公開したいと思います。とは言っても全体に公開するのではなく、知り合いだけの公開なんですが…

なぜかといいますと、今回の作品は一部のシーンで人物を使っているため、肖像権の関係で全体に公開するのはちょっとまずいんですよね。というわけでその方に許可をとって、自分の知り合いだけに公開することで了承していただいたという次第であります。ダウンロード方法は以下。だいたい一週間くらいで削除する予定です。

1. こちら(アップローダーを使用しております)のリンクからダウンロードページにアクセスします。
2. 「キーワード」の欄に自分の本名のうち、名前の部分をローマ字で入力してください。
   (ex.「佐藤 広樹」だったら"hiroki"と入力)
3. 後は画面の指示にしたがってダウンロードしてください。

なお、動画の視聴には"Quicktime"がインストールされている必要があります。"iTunes"が入っている方は自動的にインストールされるので大丈夫です。

追伸

こんなところで言うのも申し訳ないのですが、撮影に協力していただき、誠にありがとうございました。

Picasaのウェブアルバムに写真を自動的に載せるpythonスクリプトを組んでみた。

これはホント前々から作っておきたかったスクリプトの一つです。

Picasaのウェブアルバムは無料ながらにして容量が1GBととても大容量なので、クローズドに写真をアップロードする場所としてとても重宝してきました。だけどどうしても納得できないことが一つだけあって、それが写真のアップロード方法。

なんで一度に5枚しかアップロードできんのじゃい

普通アルバムに画像をアップロードしたいと思ったら2〜30枚くらいは軽くこえるだろうーどうして5枚だけ、しかもいちいち選択し直さなきゃならないのかとっても謎。

いや、「そんなことしないでPicasa使えば楽にアップロードできるのに」って思う人がいるかもしれません。そりゃそうなんですけど、個人的にあのアプリケーションは好きになれないのです。べつにキャッシュしなくてもいい画像ファイルもキャッシュしてしまうし。そんくらい自分でちゃんと整理します。

というわけで不自由なら自分でなんとかするの精神で、Picasaのウェブアルバムに指定したフォルダーの写真をすべて自動的にアップロードするpythonスクリプトを作りましたので公開します。
まず、使用するには以下の環境が必要です。

• Python 2.6
• gdata-python-client

ソースみれば大体何やってるかわかると思いますが、画像が大量に入っているフォルダを用意しておいて、
$ python picasa.py new_album ./photos/
とか指定してやれば"new_album"っていうアルバムを新しく作って、./photos/ フォルダの画像すべてを自動的にアップロードしてくれます。すごく楽ちん。

#!/usr/bin/env python
#-*- coding:utf-8 -*-

import sys, os
import gdata.photos.service

USER_EMAIL = "ログインに用いるeメールアドレス"
USER_PASSWORD = "パスワード"

def main():
 # ログイン開始
 gd_client = gdata.photos.service.PhotosService()
 gd_client.email = USER_EMAIL
 gd_client.password = USER_PASSWORD
 gd_client.source = "Picasa-AutoUploadApp"
 print u"Picasaへのログインを開始します。"
 gd_client.ProgrammaticLogin()
 
 # アルバムを新しく追加
 target_album = gd_client.InsertAlbum(title=sys.argv[1], summary='')
 print u"新しいアルバム \"%s\" を追加しました。" % sys.argv[1]
 
 # 指定されたディレクトリから写真の読み込み
 filepath = []
 for root, dirs, files in os.walk(sys.argv[2]):
  for f in files:
   if ( f.endswith(".JPG") or f.endswith(".jpg") ):
    filepath.append( (os.path.join(root, f), f) )
 
 # 写真の追加
 print u"写真の読み込み完了。アップロードを開始します。"
 
 album_url = "/data/feed/api/user/default/albumid/%s" % target_album.gphoto_id.text
 for ( i , (photo_path, photo_name) ) in enumerate(filepath):
  print u"(%i/%i) \"%s\"をアップロードしています..." % (i+1, len(filepath), photo_name),
  gd_client.InsertPhotoSimple(album_url, photo_name, '', photo_path, content_type="image/jpeg")
  print u"\t完了"
 
 print u"すべての写真のアップロードが正常に行われました。"
 
if __name__ == '__main__':
 sys.exit(main())

これでやっとひたすらボタンとにらめっこする重荷から開放される…

追記

これを書いている途中で気になったんですが、pythonってprint文で自動的に改行させないようにしたい場合、
print "hogehoge",
ってコンマを最後に付けるんですね。これなんとかなんないでしょうか？python 3.0ではさすがに修正されているとは思いますが、すんごくsyntax errorっぽい…

Hatena::Counting

「自分の大学生活は、一体あと何日で終わるのか？」

久しぶりの投稿です。なぜあまりブログを書かなくなったかというと特に意味はないんですが、強いて言うならいろいろ忙しかったため、あまりブログのネタになるようなことをしてなかったことくらいでしょうか。

はてなから新しく、Hatena::Countingというサービスが始まったようです。で、どういうサービスかというと、自分の記念日とか目標となる日から現在までのカウントダウンやカウントアップを簡単に作れるサービス。あんまりインパクトはないんですが、痒いところに手が届くようなサービスですね。

というわけで、自分もちょっくら大学生活のカウントダウンを表示して「これくらいしか時間はないんだ」と常に言い聞かせようかなと思ったんですけど、作ってみたらまだこれくらいもあるんだとか思ってしまったのは秘密です。そのリンクはこちらから。自分と同学年の方とかはちょっくら参考にしてみてください。

ヘッドホンアンプを自作してみた。

お久しぶりです。約一ヶ月ぶりの更新となりますが、べつに飽きているわけじゃなくこれからもぼちぼち続けていきます。

さて、最近友人と話していたら、どうやらその人が自作のオーディオ関係に詳しいらしく、自作のヘッドホンアンプを作っているようなんです。そこで自分もちょっくらやってみるか！ということなんで、友人の手助けをむっちゃ受けながら自分もヘッドホンアンプを自作してみました。いやぁほんとどうもありがとうございますって感じです。あとでちゃんとお礼を言わないと…

回路図はまんま友人が作ったのを流用したんでどっちかっていうと自作じゃないかもしれないですが、実際に音が出てくるとそれなりに感動です。製作費用は部品代だけで考えると2,000円くらい。市販での実売価格が7,000〜40,000円くらいなのを考えるとそれなりに良かったかなと思います。

(それにしても電子工作関係は未知の領域だったんでホント大変でした。工学部あたりだと当たり前だろとか思う知識が全くないんで手探り状態。コンデンサや抵抗の種類がこんなにあるとは思わんかったというレベルです。)

エッシャーくん、広がる。

• AS3でもエッシャーっぽい絵を出してみた
• エッシャー
• エッシャーたん
• エッシャーっぽい絵を生成する "escherize" を OpenCV で作ってみた

なんか様々な方が様々な実装方法でエッシャーくんを作っているようですね。ここまで話が広がるとは思っていませんでした！

アニメーションでぐるぐるさせている方もいまして、自分としてはとっても大満足なんですけど、やっぱりリアルタイム処理はできなさそうですね…アルゴリズム上、どうしてもサイクル数がバカ高いexp,log,sin,cosなどの関数を全ピクセルに対して行う必要がありますから、全体の計算量が半端ないことになってますし。

仮にこのコードをC/C++にrewriteして実用化させようとかいうときには、やっぱり速度とかに気をつけなければならないんでしょうね。やっぱりpythonだけでなく、もうそろそろC++について本格的に学ばないと。
(それはそうとして、gccやLLVMの最適化処理って現状ではどうなっているんでしょうね？もはや人様が下手に最適化するよりも効率が良くなったりしているんだろうか？)

なんですとー

なんですとー

いつの間にか知らないところで自分のブログのアクセス数が一気に8〜9倍近く上昇していたようで。
こんなちまちまやってる小規模ブログに一体なぜ…と思って調べてみたら、どうやら前回のエッシャーくんがオレンジニュースに掲載されたのと、はてブ経由で集まってきたのが原因のようですね。

まぁこれが過ぎればまた見る人だけ見るちっちゃいブログに戻りますし、これからもちまちま自分がおもしろいと思えるものを作っていきます。すべて世は事もなし。

追記

http://atmarkjojo.org/archives/2008/2008-09-25-001740.html

ついさっき知りました。ゴゴゴゴゴ…スタンドも月までブッ飛ぶこの衝撃…よりによってそこかぁ…

エッシャーっぽい絵を生成する「エッシャーくん」を作ってみた。

エッシャーくん

エッシャーっていう画家は知っていますか？分かんない人のために説明しますと、こんな感じのふしぎーな絵を書いている人です。名前は知らなくても一度は見たことがあるのではないでしょうか。

それでなんですが、適当な画像からなんかエッシャーっぽいふしぎな画像を生成するフィルタ「エッシャーくん（仮称）」をPython Imaging Library(PIL)で作ってみました。これを使えばどんな画像もエッシャーっぽい世界にご招待です。ソースは近々公開します。

追記(09/24)

ソースコードをアップロードしました。subversionで管理されてますので、
svn checkout http://svn.coderepos.org/share/lang/python/escher Somewhere
でチェックアウトしてください。

たとえば、こんな感じのイラストにエッシャーくんを適用させてみると…

こんな感じの絵が出来上がりました。これはどこまでズームさせても終わりがない、無限に続いている絵になっています。

この、なんかごっついリングもエッシャーくんの手によれば…
(Photo by the justified sinner, "Roofing Nail Ring 1")

どこまでも続く絵に早変わり！なんかうにょんうにょんしています。

この画像はどーなるんだろ…

なんかコーヒーっぽい…

どーやって作っているのか

エッシャーくんは"A logarithmic image transformation"っていう論文を元に作られているので、気になる人はこれを見れば大体のことが分かると思いますが、一応大まかな解説をしてみると、複素平面をlog(z)で変形して、絵が循環するように座標をスケール、回転させてからexp(z)で繋げているようです…

あーやっぱ説明が下手くそなので、見てもらったほうが早いです。

余談

本当はprocessingで作る予定だったんですけど、あまりにもprocessingの数学ライブラリがアレだったんで途中からpythonで作り直しました。(対数関数すらないのにはビックリした。そりゃああんま使わないかもしれないけどさー）

それにしてもpythonの開発速度の早さは半端ない。まさか100行未満で出来るとは思わんかったです。

リーマン破産法申請したんですって

三連休明けの東京市場やアジアの株式市場は、急落が避けられないだろう。バンク・オブ・アメリカ（バンカメ）(BAC.N: 株価, 企業情報, レポート)によるメリルリンチ(MER.N: 株価, 企業情報, レポート)買収の話しが１５日の米株市場でどの程度好感されるかわからないが、日経平均は年初来安値をトライする展開となる可能性が高いとみている。
リーマン破産法申請：識者はこうみる | Reuters

いや、たぶん全員知っているとは思うんですが、自分は前回のとは別の旅行に4日間くらい滞在して今帰ってきましたので、さっき始めてこのニュースを目にしました。やばい時代の流れに完璧に取り残されてる。
あんだけ救済するだのなんだのやっといてこの結果ですか！こりゃあさらなる信用不安が引き起こされて、株式市場は大荒れですよ。

頼みの綱であったBRICsも成長が軟化していますし、今や世界経済は下落するリスクはあれど、上昇するリスクは殆どありません。
最悪のパターンは信用不安の連鎖で「第2のリーマン」となる金融機関が次々と登場することですね。流石にそうなる前に政府が何らかの対応策を講じるとは思いますが、はたしてそれで信用が回復するかどうか…

Flight404のソースコードを読み解く part4

パーティクルを放出させる

旅行から帰ってきましたので、ようやっと更新することができました。まだまだ飽きてはいませんよ。
前回はエミッタを放出させるところまで行ったので、今回はパーティクルを等速で放出させるところまでやっていきます。

今回のソースは以下。相変わらず_gl.pdeは変わっていないので省略。新しくparticle.pdeを加えました。

Particle3.pde(draw()の関数だけ変更されているので、その部分だけ記載。)

void draw(){
  background( 0.0 );
  perspective( PI/3.0, (float)width/(float)height, 1, 5000 );
  
  gl.glDepthMask(false);
  gl.glEnable( GL.GL_BLEND );
  gl.glBlendFunc( GL.GL_SRC_ALPHA, GL.GL_ONE );
  
  pgl.beginGL();
  emitter.exist();
  pgl.endGL();
  
  if( mousePressed ){
    emitter.addParticles(20);
  }
  
}

emitter.pde

class Emitter {
  Vec3D loc;
  Vec3D vel;
  Vec3D velToMouse;
  
  color myColor;
  
  ArrayList particles;
  
  Emitter(){
    loc = new Vec3D();
    vel = new Vec3D();
    velToMouse = new Vec3D();
    
    myColor = color( 1, 1, 1 );
    
    particles = new ArrayList();
  }
  
  void exist(){
    setVelToMouse();
    findVelocity();
    setPosition();
    
    //rendering Emitter
    gl.glEnable( GL.GL_TEXTURE_2D );  //Enable texture mapping
    iterateListExist();
    render();
    gl.glDisable( GL.GL_TEXTURE_2D );
  }
  
  void setVelToMouse(){
    velToMouse.set( mouseX - loc.x, mouseY - loc.y, 0 );
  }
  
  void findVelocity(){
    vel.interpolateToSelf( velToMouse, 0.2 );
  }
  
  void setPosition(){
    loc.addSelf( vel );
  }
  
  void iterateListExist(){
    pgl.bindTexture( particleImg );
    
    for( Iterator it = particles.iterator(); it.hasNext(); ){
      Particle p = (Particle)it.next();
      if(!p.ISDEAD){
        p.exist();
      } else {
        it.remove();
      }
    }
  }
  
  void render(){
    pgl.bindTexture( emitterImg );
    renderImage( loc, 200, myColor, 1.0 );
  }
  
  void addParticles( int _amt ){
    for( int i=0; i<_amt; i++ ){
      particles.add( new Particle( loc, vel) );
    }
  }
}

particle.pde

class Particle {
  int len;
  Vec3D[] loc;
  Vec3D startLoc;
  Vec3D vel;
  
  float radius;
  float age;
  int lifeSpan;
  float agePer;
  boolean ISDEAD;
  
  Particle( Vec3D _loc, Vec3D _vel ){
    radius = random( 20, 60 );
    len = (int)radius;
    loc = new Vec3D[len];
    
    startLoc = new Vec3D( _loc.add( new Vec3D().randomVector().scaleSelf( random(5.0) ) ) );
    
    for( int i=0; i
      loc[i] = new Vec3D( startLoc );
    }
    
    vel = new Vec3D( _vel.scale( 0.5 ).addSelf( new Vec3D().randomVector().scaleSelf( random( 10.0 ) ) ) );
    
    age = 0;
    lifeSpan = (int)radius;
    agePer = 1.0;
  }
  
  void exist(){
    setPosition();
    render();
    setAge();
  }
  
  void setPosition(){
    for ( int i=len-1; i>0; i-- ){
      loc[i].set(loc[i-1]);
    }
    loc[0].addSelf( vel );
  }
  
  void render(){
    color c = color( agePer, agePer*0.75, 1.0 - agePer );
    renderImage( loc[0], radius*agePer, c, 1.0 );
  }
  
  void setAge(){
    age++;
    
    if( age > lifeSpan ){
      ISDEAD = true;
    }else{
      agePer = 1.0 - age/(float)lifeSpan;
    }
  }
}

emitter.pdeの中で特に注目すべきなのは、iterateListExist()の追加。これによって各パーティクルのレンダリング処理を行い、死んだパーティクルをリストから除外している。

particle.pde の中で特によくわかんなかったのが、loc[]の存在について。別に1個だけでいいじゃんと思っていたが、これは後に尻尾というか、パーティクルの軌跡を描画するために使うようだ。loc[0]に一番最新のパーティクルの位置が格納されていて、時が経つごとに順次入れ替えを行っている。

次によくわかんないというか、感覚的じゃないなと思ったのが、new Vec3D( _loc.add( new Vec3D().randomVector().scaleSelf( random(5.0) ) ) )などのベクトル処理。
randomVector() はベクトルをランダムな単位ベクトルに変更するメソッドで、いろいろと便利。scaleSelf()は文字通り、ベクトルを拡大したり縮小したりするメソッド。この返り値を_locベクトルに加算させていることは理解できるのだが、いかんせん感覚的じゃないのが困り者。vec1+vec2とか出来れば理解しやすいのにとか思うけど、javaは演算子のオーバーロードをサポートしていないから仕様がないのか。

Flight404のソースコードを読み解く part3

エミッタを表示させる

前回は前準備の段階で終わったので、次はパーティクルを放出させるエミッタを表示してみる。エミッタはマウスに追随するようにばね運動をする。_gl.pdeは前回と変わっていないので省略する。コードは以下。

Particle2.pde

import toxi.geom.*;
import processing.opengl.*;
import javax.media.opengl.*;

// Settings about OPENGL
PGraphicsOpenGL pgl;
GL gl;

// Settings about Emitter
Emitter emitter;

// Settings about images
PImage particleImg;
PImage emitterImg;

void setup(){
  size( 600, 600, OPENGL );
  colorMode( RGB, 1.0);
  
  hint( ENABLE_OPENGL_4X_SMOOTH );
  
  pgl = (PGraphicsOpenGL) g;
  gl = pgl.gl;
  gl.setSwapInterval(1);
  
  initGL();
  
  particleImg = loadImage( "particle.png" );
  emitterImg = loadImage( "emitter.png");
  
  emitter = new Emitter();
}

void draw(){
  background( 0.0 );
  perspective( PI/3.0, (float)width/(float)height, 1, 5000 );
  
  gl.glDepthMask(false);
  gl.glEnable( GL.GL_BLEND );
  gl.glBlendFunc( GL.GL_SRC_ALPHA, GL.GL_ONE );
  
  pgl.beginGL();
  emitter.exist();
  pgl.endGL();
}

emitter.pde

class Emitter {
  Vec3D loc;
  Vec3D vel;
  Vec3D velToMouse;
  
  color myColor;
  
  Emitter(){
    loc = new Vec3D();
    vel = new Vec3D();
    velToMouse = new Vec3D();
    
    myColor = color( 1, 1, 1 );
  }
  
  void exist(){
    setVelToMouse();
    findVelocity();
    setPosition();
    
    //rendering Emitter
    gl.glEnable( GL.GL_TEXTURE_2D );  //Enable texture mapping
    render();
    gl.glDisable( GL.GL_TEXTURE_2D );
  }
  
  void setVelToMouse(){
    velToMouse.set( mouseX - loc.x, mouseY - loc.y, 0 );
  }
  
  void findVelocity(){
    vel.interpolateToSelf( velToMouse, 0.2 );
  }
  
  void setPosition(){
    loc.addSelf( vel );
  }
  
  void render(){
    pgl.bindTexture( emitterImg );
    renderImage( loc, 200, myColor, 1.0 );
  }
}

Particle2.pdeについては、Emitterクラスと画像の初期化、並びにemitter.exist()メソッドを付け加えただけ。主な処理はemitter.pdeで行っている。

実際のレンダリングはrender()が担当している。まずbindTexture()で貼り付けるテクスチャを指定して、_gl.pdeで定義したrenderImage()で描画している。その他のメソッドではエミッタの位置をマウスに追随させるようにしている。

interpolateToSelf()はtoxi.geom.Vec3Dで定義されている、自身のベクトルを他のベクトルに置き換えるメソッド。要はvelToMouseを0.2倍したベクトルに置き換えているだけ。この処理によって、エミッタはマウスの周辺でばね運動のような動きをする。

おーし次からが本番だ。

Flight404のソースコードを読み解く part2

何も表示されないアプリケーションを作る

Flight404のサンプルではOpenGLを使用して平面に別途用意したパーティクルの画像を貼り付けて、表示するという手法をとっている。そのためP3Dと比べて、非常に高速で高精度な3D画像を描画できるのだが、その分生のOpenGLを叩かないといけないため、コードはいささか面倒になってくる。そこでとりあえずソースコードの中から何も表示されないアプリケーションを作って、コードの流れについて調べてみることにした。コードは以下。

particle1.pde

import toxi.geom.*;
import processing.opengl.*;
import javax.media.opengl.*;

// Settings about OPENGL
PGraphicsOpenGL pgl;
GL gl;

// Settings about some functions

void setup(){
  size( 600, 600, OPENGL );
  colorMode( RGB, 1.0);

  hint( ENABLE_OPENGL_4X_SMOOTH );

  pgl = (PGraphicsOpenGL) g;
  gl = pgl.gl;
  gl.setSwapInterval(1);

  initGL();
}

void draw(){
  background( 0.0 );
  perspective( PI/3.0, (float)width/(float)height, 1, 5000 );

  gl.glDepthMask(false);
  gl.glEnable( GL.GL_BLEND );
  gl.glBlendFunc( GL.GL_SRC_ALPHA, GL.GL_ONE );

  pgl.beginGL();
  // write something...
  pgl.endGL();
}

_gl.pde

int squareList;

void initGL(){
  pgl.beginGL();
  squareList = gl.glGenLists(1);
  gl.glNewList(squareList, GL.GL_COMPILE);
  gl.glBegin(GL.GL_POLYGON);
  gl.glTexCoord2f(0, 0);  gl.glVertex2f(-0.5, -0.5);
  gl.glTexCoord2f(1, 0);  gl.glVertex2f( 0.5, -0.5);
  gl.glTexCoord2f(1, 1);  gl.glVertex2f( 0.5,  0.5);
  gl.glTexCoord2f(0, 1);  gl.glVertex2f(-0.5,  0.5);
  gl.glEnd();
  gl.glEndList();
  pgl.endGL();
}

void renderImage( Vec3D _loc, float _diam, color _col, float _alpha){
  gl.glPushMatrix();
  gl.glTranslatef( _loc.x, _loc.y, _loc.z );
  gl.glScalef( _diam, _diam, _diam );
  gl.glColor4f( red(_col), green(_col), blue(_col), _alpha );
  gl.glCallList( squareList );
  gl.glPopMatrix();
}

まずhint( ENABLE_OPENGL_4X_SMOOTH )を指定して、4xのオーバーサンプリングを許可。setSwapInterval()については、主に画面がちらつく現象を防ぐ目的のようだ。その後initGL()を呼び出して、_gl.pdeで書いたOpenGL周辺のセットアップを開始する。

initGL()では、正方形のポリゴンを作って画像を貼り付けるという作業が何回も繰り返されるので、ディスプレイリストの機能を利用する。テクスチャのUV座標と頂点をそれぞれglTexCoord2fとgl.glVertex2fで指定して、正方形を描画している。

renderImage()で実際に画像の描画を行う。glPushMatrix()とglPopMatrix()で行列スタックを取得、開放し、その間に平行移動やスケールをもってくることで、ローカル座標が狂わないようにしているようだ。

次にdraw()の項を見てみる。glDepthMask(false)で陰線消去を行わないようにして、glBlendFunc( GL.GL_SRC_ALPHA, GL.GL_ONE )で画像の加算処理を指定している。そしてbeginGL()とendGL()の間にOpenGLの描画命令を組み込むことで、実際にOpenGLで描画処理が行われるようになるという次第。

Flight404のソースコードを読み解く part1

概要

Processing(Proce55ing)界で知らない人は皆無と言われる(と勝手に思っている)Flight404のRobert Hodgin。とりあえず彼のことについて説明しますと、こんな感じの映像をProcessingで作ってしまうような人です↓

Weird Fishes: Arpeggi from flight404 on Vimeo.
(悔しいことに、ルックスもイケメンだ！)

で、実は彼はブログで「簡単」な(彼にとっては簡単なだけで他の人にとっては難しいです)Proce55ingのソースコードを公開しています。が、ワールドワイドではともかく、日本語のブログで彼のコードについて言及しているサイトやブログが全くない。それで自分の丁度いい勉強にもなりますんで、適当にFlight404のSRC: Particle Emitterで気づいたことを書き留めておくことにします。

もしかしたら飽きて、中途半端なところでやめるかもしれませんが、とりあえずそんときはそんときで。てきとーにやっていきます。

Processing(Proce55ing)と流体力学

タイトルほどそんな大層なことはしていないんですが、Proce55ingと流体力学を使用して、川のような流れの中にまるーい物体を入れたときの水流の動きをアートっぽく可視化してみました。結構綺麗な模様ができて、自分としては満足です。

左から右に水が流れているもんだと思ってください。こんな感じに水流は円柱を避けて進んでいきます。

今度は円柱が時計回りに回転して、渦が発生しているようなときの水流の流れ。この渦の流れが強くなっていくと…

こんな感じにどんどんねじ曲がっていって…

最終的にはこんなかんじの、歪んだ水流になってしまいました。
本当は上の画像、アニメーションしてもっと面白い感じの映像になるんですけど、動画のキャプチャの仕方が分からないのでとりあえず画像だけ。

コードはhTakaさんのStreamDrawingの速度場を流体力学のに書き換えただけですし、しかもめっちゃ汚いのでとりあえず非公開です。興味のある方はそちらのコードを参照してください。hTakaさん、どうも有難う御座います。

それにしてもProce55ingはぱぱぱっといろんなことが出来て面白いです。なんかこれ使ってもっと有り得ないような映像できないかな。

ほんとどうでもいいですけど

(Photo by Jamoker, "The Motion of the Ocean")

ほんとどうでもいいですけどオイラーの運動方程式

に居座っている↑こいつがけっこうかわいいと思っている人はぜったい自分だけじゃないはず。
そんだけです。はいどーでもいー早く仕事しよう。

追記
と思ってたらやっぱ結構いた…
http://www.google.co.jp/search?hl=ja&q=%E3%82%AA%E3%82%A4%E3%83%A9%E3%83%BC%E3%80%80%E6%96%B9%E7%A8%8B%E5%BC%8F%E3%80%80%E3%81%8B%E3%82%8F%E3%81%84%E3%81%84

どうやらDjangoの登録方法が変更になったみたい

どうやらDjangoはバージョン0.97から管理画面にModelを登録する方法が変更されているみたいで。今までは対象のModelクラスにAdminクラスメソッドを作る方法を採用していたんですが、Modelクラスに管理画面の設定が混在せざるを得ない仕様はやっぱり完璧主義者には馴染めないものだったんですね。

で、具体的にどうするかっていうと対象のアプリケーションディレクトリに新しくadmin.pyファイルを作成して、以下のように記述します。

#!/usr/bin/env python
#-*- coding:utf-8 -*-
from blogs.models import Entry,Blog,Comment
from django.contrib import admin

##
## The Configuration for Admin.
##
class EntryInline(admin.StackedInline):
	model = Entry
	extra = 1

class BlogAdmin(admin.ModelAdmin):
	inlines = [EntryInline]
	list_display = ('title','pub_date','user')
	list_filter = ['pub_date']

# registration
admin.site.register(Blog, BlogAdmin)

admin.site.register()の部分で管理画面の登録を行ない、具体的な設定は専用のクラスを新しく作成することで解決します。
しかしこのままではadmin.pyが呼び出されないので、プロジェクトのurls.pyの先頭に以下の2行を加えます。

from django.contrib import admin
admin.autodiscover()

これでadmin.pyが呼び出され、無事管理画面の登録が行なわれるようになりました。

Djangoのオンラインドキュメントではmodels.pyに直接Adminの記述を行なっているんですが、これだとviews.pyとかがmodels.pyをインポートしたときに再びadmin.site.register()が実行されてしまいますので、二重登録のエラーが生じてしまいます。モデルと管理画面の設定を分離するという意味でも、こちらのほうがより便利じゃないかなーと思います。

Django始めました。

っていうとなんか冷やし中華っぽいですよね。

というわけで、新しくDjangoの勉強をすることにしました。ここでジャズ好きの人は「ラインハルトの勉強って…ジャズでも始めたんかい」って思うかもしれませんが、このDjangoはpythonで動くWebアプリケーションのフレームワークなんです。

Django（ジャンゴ）は、Pythonで実装されたWebアプリケーションフレームワーク。 model-view-controller デザインパターンを緩やかに従う。 もともとはローレンス (カンザス州)にある World Companyのために、ニュース系のサイトを管理する目的で開発され、2005年7月に BSD License で公式にリリースされた。フレームワークはバンド gypsy jazz のギタリスト ジャンゴ・ラインハルト にちなんで命名された。(From Wikipedia)

あぁなんてマニアックなんでしょう。すいませんこんなつまんないブログで。でもやってみたかったんです。とりあえず勉強してるんですから、なんか成果を出すことができたらいいですね。

追記
あーでもジャズもなんか面白そうですよね。やってみたい。

ETF投資を勧めるのはいいけど、数学的ってつけるのはちょっと勘弁

なんか最近ETFと国債のポートフォリオ投資がとても人気なようで、あちこちにETF投資を勧める記事が散見されます。いちいち数千もある資産の変動を考える必要もないですし、なによりCAPM理論によればETFなどのインデックスファンドと国債の組み合わせが最も合理的でパフォーマンスのよいポートフォリオを組めると数学的に証明されているので、そのことも人気の一つとしてあるでしょう。

うーん…別にどうというわけではないんですけど、ETFを勧めている人はもしかしてCAPM理論を未だに信じているんでしょうか？CAPMを使って現在価値の計算をやっている人はもう殆どいないというこのご時世に…
というわけで、とりあえず自分の考えを整理するという意味でも、CAPM理論について書いてみることにしました。

大前提

ETFと国債の組み合わせによるポートフォリオが「数学的に合理的」なポートフォリオとなるためには、以下の4つの仮定を満たすときに限られます。

1. 全世界のすべての投資家が『マーコヴィッツの平均-分散モデル(mean-variance model)』を元に投資を行なっている。
2. 投資家のすべてがすべての投資可能な資産について、その期待収益率、分散、共分散に同じ値を付与している。
3. 無リスク資産の貸し借りに関する一意的な無リスク資産が存在し、取引コストは存在しない。
4. 投資期間は1期間であるとする。

それではこれらの仮定をもう少し深く見てみることにしましょう。

1.全世界のすべての投資家がマーコヴィッツの平均-分散モデルを元に投資を行なっている。

『マーコヴィッツの平均-分散モデル』というのは全てのポートフォリオの組み合わせの中から最もリスクが低くなるポートフォリオの組み合わせを選択するモデルです。確かにこのモデルが主張していることは、リスクを最大限回避するという意味で重要かもしれません。ですがこのモデルを全ての投資家が利用していると仮定するのには流石に無理があると思います。というか個人投資家でETFは知っていても、このモデルを知っている人は少ないですよね。主従逆転状態です。

2.投資家のすべてがすべての投資可能な資産について、その期待収益率、分散、共分散に同じ値を付与している。

結構きつい仮定です。分散、共分散については過去の値動きから算出できるのである程度同じ値になる可能性はありますが、期待収益率については過去の値動きがまったく当てにならないので、ファンダメンタルズから推測するしかありませんし、その推測も個人でバラバラです。それなのに投資家の全てが同じ期待収益率を推測していると仮定するのはかなりナンセンスな話です。

3.無リスク資産の貸し借りに関する一意的な無リスク資産が存在し、取引コストは存在しない。

これはある程度納得できる仮定です。国債などの流動性が高い無リスク資産においては買値と売値の幅は非常に狭いですし、取引コストも扱う金額が大きければ微小量として無視できるくらいの量です。

4.投資期間は1期間であるとする。

結構この仮定を忘れられている気がします。ETFの投資が合理的になるときの条件は、『1期間(1年間とか1ヶ月とか)』で取引した場合の話であって、多期間で投資した場合では異なった結果が得られます。どういうことかというと、この仮定では1期間の間に新しく売買取引を行なうことができないのです。たとえ暴落、暴騰が起こってもです。理論を簡単にするためとはいえ、現実的でない仮定です。

というわけで

CAPM理論は結構無理な仮定をおいた上で証明してしまった理論なので、これらの仮定を1つでも信じられないという人はあまり使えない理論です。

ただインデックスファンドとアクティブファンドの平均的なパフォーマンスを比べてみると、インデックスファンドのほうが殆どのアクティブファンドよりも良いパフォーマンスを示したということは事実ですので、困ったらETFに投資してみるといいかもしれません。

まとめ

ETFと国債への投資はパフォーマンスは良いかもしれないけど、数学的に信用できるものではない。投資方法の一つとして見るべき。

追記(8/10)
対象の記事にコメントできませんね。なんか本人の知らないところでうだうだ言うのはあんまよくないと思っているので、引用部分を置き換えました。

原油先物価格が順調に

暴落しています。7月9日の北海道・洞爺湖サミットで、原油価格高騰の解消に向け産油国へ増産を促すことが合意されたことをきっかけに、先物の売りが相次いでいるようです。

米商品先物取引委員会(CFTC)は今回の上昇の原因を単純な需給不均衡と見ているようですが、それだけではここまで急激な一方向の上昇(1年半で2倍弱！)を説明することはできません。少なくとも原因の一つとして投機が絡んでいたことはほぼ確実だと思っていますので、まだまだ価格は下がると予想しています。早くガソリン代が安くなるといいですねぇ。

Googleのストリートビューが凄すぎる

題名の通り。いつの間にかGoogle マップにストリートビュー機能が追加されていたみたいなので試してみたら、予想以上にインパクトのある機能でした。

何がそんなに凄いのかっていうと、カバーする地域がとても広いところ。右上の「ストリートビュー」をクリックすると道路が青く光るので、見たい地点をクリックしてストリートビューを表示します。画面をドラッグして360度自由自在に動かすことも可能。

(上は自分が住んでいる仙台で最も人通りが多いアーケード街の画像。青い道路ならばどこでもストリートビューで見ることができます。)

(フォーラスが改装中ってことは、7月の中旬か下旬に撮られた？)

いやぁーそれにしてもGoogleは本当にヤバい。自分が普段通っている道路が鮮明に映されるときのインパクトは相当です。
ストリートビューで自分の住んでいるマンションが映されたときはおもわず笑ってしまいました。

(ついでにパリの凱旋門のストリートビューを紹介します。とてもきれい。)

linuxで一斉に画像をリサイズしたい場合

とある事情で2,300枚ある画像をすべて800x600のjpgに変換しなきゃなんなくなりました。

最初は面倒くさいから「縮小専用。」をwineで動かそうかなっていう安易な発想をしていたんですが、どうもドラッグ&ドロップのイベントが効かなくなっているようで動く気配がありません。しょうがないのでlinux上でなんとかしなきゃということで調べてみたら、mogrifyコマンドがどうやら有効なようで。対象のディレクトリに移動してから、
$ mogrify -geometry 800x600 -quality 80 *.jpg
で一発でした。ただし画像はすべて上書きされてしまうので、不安だったらバックアップをとってください。
このmogrifyコマンド、スケールだけじゃなくてどうやらクリップやドロップシャドウなんかもできるらしいので、結構高性能です。

それにしてもlinuxは覚えると圧倒的に楽ですが、それまでが大変ですね。

The ting tings - Great DJ & Shut Up And Let Me Go

Motion Graphicsとしてはあまり真新しいものではないですが、最近気に入っているバンドなので掲載しました。

Wowowをぼーっと見ていたら偶然このバンドのPVが流れていたので、気になって後で調べてみたらびっくり。このバンドってiPodのCMに使われたことがあるみたいです。The fratellisの件といいDaft Punkの件といい、iPodのCMに使われるアーティストと自分が好きなアーティストにはどっかしら共通点があるのかもしんないです。

The Ting Tings - Great DJ

これが最初にwowowで見たPV。この人たちは本当にやる気があるんでしょうかというくらいダンスが適当。けどそれがいい。曲調と遊び具合がすんごく木村カエラっぽいなぁと思ったら、やっぱりカエラ自身もそう思ったらしく、お勧めだよと言ってました。

The Ting tings- Shut up and let me go

手の中に映像があってさらにその手の中に次の映像があって…というよく使われている手法ですが、自分はこの手法が大好きです。途中のカンフーはたぶんアクションシーンが好きで入れたんだと思います。だって唐突でよくわかんないし。

ちなみに名前の由来は中国語のTing(聴く)からきているらしいので、そういう変な意味はないです。

CodeReposをやることにしました。

CodeReposのcommit権作っといたよーっいうメールが届いてました！

いままでプログラムを書く際にはローカルなレポジトリを作ってそこで管理するようにしていたんですが、これだと学校などで思いついたアイデアを実装したいときちょっとめんどいですし、なによりネット上で常に最新のリビジョンを公開することができません。
それでいろいろと調べてみた結果、CodeReposが一番自分の目的に合っているようなサービスだと感じましたので、これからはこのサービスを使って開発をしていきたいと思います。

それにしても、このサービスは基本的にハイレベルなプログラマーがちょこちょこcommitしているような印象がありますので、ちょっと緊張しますね。

Black-Scholes Equation in Python 改良版

前回のコードをある程度改良しましたので、載せておきたいと思います。主な改良点は以下。

• ヨーロピアン(コールorプット)オプションのデルタ、ガンマ、ベガ、シータ、ローを求める関数を追加しました。
• ヨーロピアン(コールorプット)オプションのインプライド・ボラティリティをニュートン法より求める関数を追加しました。
• その他若干の名称変更や、引数に関する説明を加えました。

数式は"Black-Scholes (1973) Option Pricing Formula"を参考にしました。

インプライド・ボラティリティに関しては適切な引数を加えないと、負の値や無限大に発散することがあるので注意が必要です。とはいえ実際の市場価格と理論価格がそこまで乖離することもないでしょうから、大抵正常な値を返すはずです。

商用、非商用に関わらず自由に使用、改変していいですし連絡する必要もないですが、無保証です。使う人はあんまいないと思いますが、一応建前として。

#!/usr/bin/env python
#-*- coding:utf-8 -*-
"""
Functions for calculating the Black-Scholes equation.
"""
import math

A1 = 0.319381530
A2 = -0.35653782
A3 = 1.781477937
A4 = -1.821255978
A5 = 1.330274429
RSQRT2PI = 0.3989422804
GAMMA = 0.2316419

def BS_Call(S,K,T,R,V):
	"""Calculate the Black-Scholes formula for a call option.
	
	Arguments:
		S:the price of the underlying stock
		K:the strike price
		R:the continuously compounded risk free interest rate
		T:the time in years until the expiration of the option
		V:the volatility for the underlying stock
	"""
	d1,d2 = _getD1D2(S,K,T,R,V)
	ret = S*_normalDist(d1)-K*math.exp(-R*T)*_normalDist(d2)
	return ret

def BS_Put(S,K,T,R,V):
	d1,d2 = _getD1D2(S,K,T,R,V)
	ret = -S*_normalDist(-d1)+K*math.exp(-R*T)*_normalDist(-d2)
	return ret

def BS_Call_IV(S,K,T,R,Pr,HV,maxNum=100,i=1):
	"""Calculate the Implied Volatility using Newton's method.
	
	Arguments:
		Pr:the real price of call option
		HV:the Historical Volatility
		maxNum:the number of repeating calculation
		i:index(Do not specify this parameter)
	"""
	if i > maxNum:
		return HV
	else:
		newV = HV - (BS_Call(S,K,T,R,HV)-Pr)/BS_Vega(S,K,T,R,HV)
		return BS_Call_IV(S,K,T,R,Pr,newV,maxNum,i+1)

def BS_Put_IV(S,K,T,R,Pr,HV,maxNum=10,i=1):
	if i > maxNum:
		return HV
	else:
		newV = HV - (BS_Put(S,K,T,R,HV)-Pr)/BS_Vega(S,K,T,R,HV)
		return BS_Put_IV(S,K,T,R,Pr,newV,maxNum,i+1)

def BS_Call_Delta(S,K,T,R,V):
	d1,d2 = _getD1D2(S,K,T,R,V)
	return _normalDist(d1)

def BS_Put_Delta(S,K,T,R,V):
	d1,d2 = _getD1D2(S,K,T,R,V)
	return _normalDist(d1) - 1.0

def BS_Gamma(S,K,T,R,V):
	d1,d2 = _getD1D2(S,K,T,R,V)
	return _normalDist(d1)/(S*V*math.sqrt(T))

def BS_Vega(S,K,T,R,V):
	d1,d2 = _getD1D2(S,K,T,R,V)
	return S*_normalDist(d1)*math.sqrt(T)

def BS_Call_Theta(S,K,T,R,V):
	d1,d2 = _getD1D2(S,K,T,R,V)
	return -S*_normalDist(d1)*V/(2*math.sqrt(T))-R*K*math.exp(-R*T)*_normalDist(d2)

def BS_Put_Theta(S,K,T,R,V):
	d1,d2 = _getD1D2(S,K,T,R,V)
	return -S*_normalDist(d1)*V/(2*math.sqrt(T))+R*K*math.exp(-R*T)*_normalDist(-d2)

def BS_Call_Rho(S,K,T,R,V):
	d1,d2 = _getD1D2(S,K,T,R,V)
	return K*T*math.exp(-R*T)*_normalDist(d2)

def BS_Put_Rho(S,K,T,R,V):
	d1,d2 = _getD1D2(S,K,T,R,V)
	return -K*T*math.exp(-R*T)*_normalDist(-d2)

def _getD1D2(S,K,T,R,V):
	vSqrtT = V*math.sqrt(T)
	d1 = (math.log(float(S)/K)+(R+(0.5*V*V))*T)/(vSqrtT)
	d2 = d1 - vSqrtT
	return (d1,d2)

def _normalDist(x):
	if x >= 0:
		k = 1.0/(1.0 + GAMMA*x)
		cnd = RSQRT2PI*math.exp(-0.5*x*x)*(k *(A1+ k*(A2+ k*(A3 + k*(A4 + k*A5)))))
		ret = 1.0 - cnd
		return ret
	else:
		return 1.0 - _normalDist(-x)

て、そんなことやってる場合じゃない。早くテスト対策行なわないと。

Black-Scholes Equation in Python

Photo by Martinlu, "Math"

やっぱ金融学んでるんだったら、オプションの理論価格ぐらい自分で出さないといけないよねー
ってわけで、あっちの方面では超がつくほど有名なBlack-Sholes Equotionをノリで実装してみました。特にこれといって特別なことはしてないですし、もう既にどっかの誰かがやってそうですが、こんくらい自分でやんないと。

#!/usr/bin/env python
#-*- coding:utf-8 -*-

import math

A1 = 0.319381530
A2 = -0.35653782
A3 = 1.781477937
A4 = -1.821255978
A5 = 1.330274429
RSQRT2PI = 0.3989422804
GAMMA = 0.2316419

def Black_Scholes_Call(S,K,T,R,V):
 d1,d2 = _getD1D2(S,K,T,R,V)
 ret = S*_normalDist(d1)-K*math.exp(-R*T)*_normalDist(d2)
 return ret

def Black_Scholes_Put(S,K,T,R,V):
 d1,d2 = _getD1D2(S,K,T,R,V)
 ret = -S*_normalDist(-d1)+K*math.exp(-R*T)*_normalDist(-d2)
 return ret

def _getD1D2(S,K,T,R,V):
 vSqrtT = V*math.sqrt(T)
 d1 = (math.log(float(S)/K)+(R+0.5*V*V)*T)/(vSqrtT)
 d2 = d1 - vSqrtT
 return (d1,d2)

def _normalDist(x):
 if x >= 0:
  k = 1.0/(1.0 + GAMMA*x)
  cnd = RSQRT2PI*math.exp(-0.5*x*x)*(k *(A1+ k*(A2+ k*(A3 + k*(A4 + k*A5)))))
  ret = 1.0 - cnd
  return ret
 else:
  return 1.0 - _normalDist(-x)

ちなみに累積正規分布関数のところは5次の項まで近似した関数を使ってます。

そういうこともある。

うーん…最近はなんかずっと『三歩進んで二歩下がる』を地で行っているような気がします。
『調べる→組んでみる→お、これ凄いんじゃね！？→あー間違ってた…成果微妙…』
の繰り返し。そりゃあ一朝一夕で全てを成すことができるとはとても思っていませんが、高揚の次に必ず待っている失望がかなりやんなっちゃう感じです。

http://ymatsuo.com/
最近見つけたサイト。「業績」の欄を見てみると、人工知能とWebの特性を生かした面白い研究をなされているようです。
日本語でこういった最新の論文を見ることが出来るというのはかなり珍しいですので、これから要チェック。

これは…

http://jihan.sblo.jp/
KOREHA…ここまで衝撃を受けたブログはもしかして初めてかもしれない。こういう地道な積み重ねは自分にも必要だと思う。

FXのレバレッジについて。

Photo by cobalt123, "Money, Hanging On"

レバレッジとは？

外国為替証拠金取引であるFXは、証券会社に証拠金を手渡すことによって、それを元手に何倍、何十倍もの通貨を取引することができ、その仕組みを「レバレッジ」と呼びます。

例えば1万円を証拠金として1ドル100円の通貨を1000個購入して、後に1ドル120円になったときに売却したとすると、購入したときの代金は100*1000 = 10万円。それが後に12万円になっているので、差益は12万円-10万円=2万円。レバレッジは10万円 / 1万円 = 10倍となります。要は金利のない借金みたいなものです。実際にかけられるレバレッジの量は証券会社によって違いがありますが、だいたい400倍くらいが最高のようです。

さて、他のFXについて言及しているブログではレバレッジについてあんまり詳しく書いていなくてちょっと不満だったので、これを数式にして、そのあとそれから生じるリスクについて考えてみることにしました。

レバレッジの数式

あなたは優秀な投資家で、常に安定した複利で資産を成長させることができるものします。それにレバレッジを常に一定の倍率でかけた場合、資産はどのように成長するのでしょうか？

まずN0を初めに投資する資産金額、rを1期間における複利、Lを1期間にかけるレバレッジの倍率と定義します。そうすると1期間後における資産金額N1は以下のように表せます。

これは第2項目までが1期間によって生じた収益、第3項が元の金額となっています。同様にこの数式を2期間後にも適応してみましょう。

このように綺麗な数式となりました。これを第n期間まで拡張すると、

となります。単純にレバレッジをかけずに運用した場合、L=1となるので、この式はrをL倍した複利で運用することに相当します。

で、結局なんなのさ？

具体例で考えてみることにします。もしあなたが年利10%で安定した運用を行えたとすると、これをレバレッジをかけずに10年間続けた場合の資産の倍率はL=1,r=0.10より
(1+0.1)^10 = 2.59 倍
となります。それではレバレッジを仮に3倍でかけたとすると、L=3となるので、倍率は
(1+3*0.1)^10 = 13.79 倍
となります。全然違いますね。Lの値がN0ではなく、指数関数に入っているところが恐ろしいです。

でもそんなおいしい話が転がっているわけがない

そりゃあ安定して10%も運用できればいいですが、普通は下落するリスクも併せ持っていることが普通です。14倍になるってことは資産が14倍下落するリスクも同様に併せ持っているということです。

仮にrが確率変数で、平均がµ,標準偏差がσだとすると、確率変数の性質よりLrの平均はLμ,標準偏差がLσとなってしまいます。このことよりリスクもリターンもLに比例していることがわかります。

しかもFXは一度でも資産があるデッドラインまで下がったら自動的にロスカットが実行されてしまいます。つまりレバレッジを上げるということは、資産を失う確率も高めるということです。

どれくらいのレバレッジなら許容できる？

追記。それではどれくらいのレバレッジならば、ロスカットとならずに比較的安定した運用ができるのでしょうか？

通貨をUSDJPYで固定して、1年間の運用を行うものと考えてみます。仮にrが正規分布に従う確率変数として、1年間のリターンをμ[%],1年間の標準偏差(ヒストリカルボラティリティ)をσ[%],ロスカットとなる利率を-d[%]とすると、限度であるレバレッジLは

となります。試しに計算してみましょう。スワップ派のサイトによると、USDJPYのHVはσ=6.95%であるので、今回は確実に安定して運用するために95%の確率で収束する2σを採用します。またμの値はUSDJPYのスワップポイントである、2.1%を採用します。そしてdを50%とすると、Lの値は
L = 50 / (6.95*2 - 2.1) = 4.24[倍]
となるので、レバレッジが4倍以内であれば、およそ95%の確率で安全圏に収束するということが分かります。

まとめ

• 確実に利益を出せるという自信がない限り、あまりレバレッジを極端にかけないほうがよい。
• 逆に言えば、確実に利益が出せる自信があるなら積極的にレバレッジをかけたほうがよい。
• 何も考えずに400倍とかのレバレッジでやっている人はただのギャンブラー。
  

雑感。

Photo by oskay, "WhiteMenorah.jpg"

Webプログラミングに関しては数多くのフレームワークが存在して(DJango, CakePHP, Ruby on Rails, etc...)いるのにも関わらず、システムトレードに関しては全くといっていいほどフレームワークが存在していない。存在していてもなぜかアプリケーション上で動くわけの分からない独自言語(MetaTraderとか)だったりで、非常に汎用性が低いのが現実です。なんでだろ？そっちのほうが使いやすいから？明らかに開発コストもかかるし、ライブラリも充実しない気がします。

うーん…これからはシステムトレードも独自言語じゃなく、フレームワークで管理する時代だと思います。これは今大雑把に考えている構想ですが、整理するために文章に直してみます。

バックテストやシステム評価はフレームワークが総合的に管理して、ユーザはシステムの開発それだけに専念すればいい。そしてそのシステムは手直すことなく、そのまま実際に売買を行うことができます。SimulationとRealの差異はフレームワークがある程度吸収します。

ここまではMetaTraderが実現していることですが、違うのはシステムの根幹となる言語をpythonで作ることでOSを気にすることなくマルチプラットフォームで動くようにして、なおかつpythonで培われた既存のライブラリを利用できるという点。まさに「巨人の肩の上に立つ」ことで、従来のアプリケーションには出来なかった領域を開拓していくフレームワークが出来ればなぁと考えています。

どうだろ？とりあえず自分のためにいろいろ勉強してみますが、もしうまく開発が進んでいったらオープンソースでリリースしてみたいです。いままで散々オープンソースコミュニティにはお世話になってきましたから、すこしでも恩返しができればいいですね。

Pythonのソートについて: sort()とsorted()

Photo by Fort Photo, "Fountains of Light"
自分のためにメモ。pythonでソートを行おうとする場合、リストに組み込まれている関数list.sort()と、単独で使える関数sorted()の2つがある。基本的に取る引数は同じなので機能自体はあまり変わりないのだが、以下の2点が異なる。

• sort()関数は返り値を出さずに、リストの本体自体を変更してしまう。これはreverse()関数と同じ。
• sort()関数はソートの対象がリストであるのに対し、sorted()関数は対象がイテレータとなる。

試しにpython上で確認してみる。

>>> x = [1,6,3,8,4]
>>> y = [1,6,3,8,4]
>>> y.sort()
>>> y
[1, 3, 4, 6, 8]
>>> sorted(x)
[1, 3, 4, 6, 8]

きちんと両方ともソートが行われているようだ。それでは次に引数cmp, key, reverseについて調べてみる。
"cmp"キーワードには、二つのオブジェクトxとyがあったときに、どちらが大きいのかを示す関数、あるいはλ式を指定する。普通はそのままで十分であるが、比較用のオブジェクトがちょっと特殊な場合に役立つ。
もしcmpキーワードを指定しなければならないような場合、通常ならばλ式の中にcmp(x,y)関数を導入する。cmp(x,y)関数はxとyを比較して、指定すべき適切な値を返してくれる関数。これを導入すれば、
『あれ？cmpキーワードってx < yの場合にはx - yを指定すればいいんだっけ？それともy - xだっけ？』とかいうような余計な心配をしなくて済む。
"key"キーワードには、比較対象とするオブジェクトはどのオブジェクトであるのかを示す関数、あるいはλ式を指定する。何も指定しなかった場合、配列やタプルであれば一番目(インデックスが0)のオブジェクトが比較対象であると判別されるようだ。
"reverse"キーワードはソート結果が昇順であるか降順であるかを指定する。Trueを指定すると降順になる。sort()してからreverse()などのように2つメソッドをくっつける必要はない。
それでは今までの流れを踏まえて試しに簡単なプログラムを組んでみる。辞書には{"名前":"点数"}などのように名前と点数のペアが保存されている。これを点数が高い順に並び替えるというプログラムは以下のように記述できる。

>>> dic = {"Hanako":"50", "Taroh":"80", "Ichiroh":"30", "Yoshio":"60"}
>>> sorted(dic.items(), cmp=lambda x, y:cmp(int(x), int(y)), key=lambda x:x[1], reverse=True)
[('Taroh', '80'), ('Yoshio', '60'), ('Hanako', '50'), ('Ichiroh', '30')]

ちなみに上の式は以下のようにも記述できる（というかそっちのほうが普通）。

>>> sorted(dic.items(), key=lambda x:int(x[1]), reverse=True)
[('Taroh', '80'), ('Yoshio', '60'), ('Hanako', '50'), ('Ichiroh', '30')]

6/30 追記
最初のソート部分が間違っていたので修正しました。以前の書き方だとxとyは同一の配列オブジェクトを参照していることになりますね。すいません。

GECCO

Photo by oldbones, "Evolution found in the trash."

GECCO は進化計算の分野で世界最大であり，最も伝統のある国際会議です． 口頭発表は投稿論文の50％以下に絞られますが，通常の発表論文とは別にポスターセッションでの発表， ワークショップでの発表，Late Breaking Paperとしての発表があります． また，学会参加者は，ワークショップやチュートリアルに無料で参加できます． この国際会議では，GAやGP，ESなど分野ごとにBest Paper Awardが授与されますが，その候補は事前の 審査で決定され，同一セッションで発表が行われます．
(http://www.ie.osakafu-u.ac.jp/~hisaoi/news/news_gecco.html)

GECCOという進化計算に関わる学会があるということを最近知りました。とりあえず一通りGA関連の論文を読んでみたいのですが、バックグラウンドの知識が殆どないので暫くは引用先、書籍とにらめっこしながら読むことになりそうです。

でも教授に直接聞いたほうが早いかも。近くにGAとかデータマイニング関係に詳しい人いないかな。

そうだ。海へいこう。

全然いままでの内容とは関係ないのですが、今日急に海を見たくなったので、海を見に出かけることにしました。自転車で。出発地点は仙台駅からで、走行距離は片道約12kmくらいです。

最初は国道線をずっと東に進むだけだったのでけっこう大きめの店（ドンキホーテとかパチンコとかブックオフとか地方の国道線によくある店）が立ち並んでいたのですが、途中から抜けて一面の小麦畑へ突入。天気も風もちょうど良くて気持ちがいい。

小麦畑が終わった後は森の中へ。森の中は自分以外誰も居ないため生活音が全くなく、とてもしずかで、ちょっと暗い。夜立ち入ったら普通に怖いと思います。なんか落書きがあちこちにあるし。
しばらく歩いて、もしかしたらこのルート間違ってるんじゃないの？って思ったあたりで謎の階段が。一生懸命のぼっていくと、

「そのとき突然、魚の目をした男は、自分がどこに向かって歩いているのか、理解した----はっきりと理解した。心臓が止まらんばかりに。
　海だ。」
ーミヒャエル・エンデ　『鏡の中の鏡』ー

いやーほんと人並みな感想で申し訳ないんですが、ほんとうに大きいです。でかい。そりゃ今まで海を見たことがないわけじゃないですが、今回はなんか感じるものがありました。うん。すごかった。
予想以上に綺麗だったので自分以外に2,3人くらいしかいない宮城の砂浜の上で、しばらく海の音を聞きながら30分くらいぼーっとしていました。

その後はしばらく自転車でうろうろ。なぜか途中のサイクリングロードでおばさんと談笑したり、帰りの国道沿いの店によってみたりとか。

最後は偶然ドスパラ（仙台にあるとは思わなかった）を見つけて、立ち寄ったらメモリが1Gで2000円という安めの価格だったので思わず購入して、そのまま帰宅。思ってた以上に距離があってけっこう疲れたり、結局なにしたかったのかよくわかんない一日でしたが、なんだかんだで楽しかったんで結果オーライです。

大きな地図で見る

MeCab+pythonという組み合わせ

photo by s1m0ne "Ocean tender love cocktail"

最近は日本語を解析してなんかできないかなぁと思っているので、形態素解析用のソフトウェアMeCabをインストールして、いろいろ遊んでいる状態です。インストールはそんなに難しくないと思いますが、とりあえず載せておきたいと思います。

まずはMeCabとその辞書、ライブラリ一式のインストール。ソースからビルドしてインストールするのでもいいですが、今回は手軽にapt-getを使用。
$ sudo apt-get install mecab mecab-ipadic mecab-utils libmecab1 libmecab-dev
で簡単にインストールできました。
ですがこのままでは辞書の文字コードがeuc-jpであるため、utf-8が標準のubuntuでは少し使いづらい。調べてみたら、IPA辞書をeuc-jpからutf-8に変えるコマンド
$ sudo /usr/lib/mecab/mecab-dict-index -d /usr/share/mecab/dic/ipadic
-o /var/lib/mecab/dic/ipadic -f euc-jp -t utf-8 -p
がありました(From ubulog)ので、それを実行。ためしにコマンドラインから

$ mecab
隣の客はよく柿食う客だ
隣 名詞,一般,*,*,*,*,隣,トナリ,トナリ
の 助詞,連体化,*,*,*,*,の,ノ,ノ
客 名詞,一般,*,*,*,*,客,キャク,キャク
は 助詞,係助詞,*,*,*,*,は,ハ,ワ
よく 副詞,一般,*,*,*,*,よく,ヨク,ヨク
柿 名詞,一般,*,*,*,*,柿,カキ,カキ
食う 動詞,自立,*,*,五段・ワ行促音便,基本形,食う,クウ,クウ
客 名詞,一般,*,*,*,*,客,キャク,キャク
だ 助動詞,*,*,*,特殊・ダ,基本形,だ,ダ,ダ
EOS

無事に文字化けせずに形態素解析が行われているようです。面白い！

次にpython上からmecabを動かしたいので、mecab-pythonバインディングをインストールします。
sourceforge.jpからmecab-python moduleをダウンロードして、メインディレクトリ上から
$ python setup.py build
$ sudo python setup.py install
でインストール。ためしにpython上から、
$ python
>>> import MeCab
>>> m = MeCab.Tagger("-Ochasen")
>>> print m.parse("羊をめぐる冒険")
羊 ヒツジ 羊 名詞-一般
を ヲ を 助詞-格助詞-一般
めぐる メグル めぐる 動詞-自立 五段・ラ行 基本形
冒険 ボウケン 冒険 名詞-サ変接続
EOS

無事に出力されていることがわかります。

あとは文章のキーワードを抽出するようなプログラムを自分で書くだけ…と思っていたのですが、なんともう実際にモジュールとして提供しているサイトを発見。しかもApache License v2.0なのでとても使いやすい。マジっすか！ということでこれも試用してみます。

まずはhttp://tanashi.s240.xrea.com/mword/にアクセス。mword.pyとして保存。ただしこれは辞書がeuc-jpであることを前提としているので、エンコーディングをutf-8にした上で、"coding: euc_jp"から"coding: utf-8"に変更します。そしてsample.pyを手本にしながら、試しにあるニュース記事のキーワード抽出を行ってみると、

キーワードを抽出させる元の文章：
「４月半ばからレンジ相場が続いてきたドル／円<JPY=>が１０日、テクニカル上の上値抵抗線だった１０６円半ばを上抜け、３カ月ぶり高値を更新した。バーナンキ米連邦準備理事会（ＦＲＢ）議長の発言などをきっかけに、市場では米当局が信用リスク問題をにらんだ金融緩和から、インフレとドル安警戒に軸足を移したとの見方が強まっている。ドルは短期的に２月高値の１０８円台への上昇を見込む声も上がっている。....」

$ python test.py | head -n 10
ドル 63.639610
市場 17.320508
通貨 16.000000
インフレ 14.000000
上昇 13.856406
姿勢 12.000000
発言 10.392305
米当局 10.182920
豪ドル 9.797959
ドル安 8.239069

と、どうやらドル通貨に関する記事であるということがわかります。きちんとキーワード抽出を行えているようです。

しかし、殆どソースコードを書かずに日本語解析まで出来てしまうと、逆に心配になってきます。
『車輪の再発明は極力行わない』というのがプログラミングでの大原則ですが、ここまで簡略化されてしまうとなんか申し訳ないですね。開発者に感謝。

Monster Television Commercial: Stork

ちなみにMonster.comっていうのはアメリカの求人サイトなので、そこらへんを考慮に入れれば何を伝えたいのか分かりやすいんじゃないかと思います。それにしてもCMとは思えないくらいの綺麗な映像です。製作はFramestore。

ついでにもう一つ。

ランダムフォレスト(Random Forest)法を用いて株価予測を行ってみた。-03

photo by Desktopography

第三回。それでは今回行った手法について概要を説明していきます。

概要

まず株価についてのデータが欲しかったので、pythonで株価をネット上からダウンロードして結合するスクリプトを作成。そしてその時系列データを元に前日のデルタ、ヒゲの長さ、出来高などの実測値、Moving Average,MACD,RSI,Bolinger Band,RW%R...などのテクニカル指標を合わせて70個ほど算出し、CSVデータとして保存します。

CSVデータをOpenOffice.orgで開いた状態。

そのデータセットを元にR上で学習させます。予測させる数値は『翌日の終値-翌日の始値』とします。一応Rのソースを下に載せておきます。

> library("randomForest") #RFパッケージの読み込み
randomForest 4.5-25
Type rfNews() to see new features/changes/bug fixes.
> d <- read.csv("dataset/2502_R.csv", header=T) #CSVデータの読み込み
> dim(d) #次元の確認
[1] 1794   83
> d.learn <- d[1:1500,] #データを学習用と
> d.test <- d[1501:1794,] #テスト用に分割する
> d.rf <- randomForest(Result~., data=d.learn, ntree=300) #RFで学習させる。
> d.pred <- predict(d.rf, newdata=d.learn) #予測。試しに学習用のデータを入力してみる
> cor(d.pred, d.learn$Result) #相関係数を測定
[1] 0.9693442 #正しく分類が行われていることが分かる
> d.pred <- predict(d.rf, newdata=d.test) #テスト用のデータを用いて予測
> cor(d.pred, d.test$Result)
[1] 0.0266478 #ほぼ無相関の値を示す

結果

どうやらきちんと学習が行われているようなのに、未来のデータに対して予測を行わせるとほとんど無相関、つまりランダムの数値を返してしまうようです。この後も目的変数を『5日後の終値-翌日の終値』としたり、さまざまな銘柄に対して学習と予測を行ってみましたが、良くても相関係数が.2くらいで、実際に使えそうなデータはありませんでした。これは「過去の時系列データを元にして株価がどれくらい動くか予測することはかなり困難である」ことを表しています。

ウィーク型の効率的市場仮説によれば、時系列データを元に予測を行うことは不可能であると言っているわけですし*1、このような結果が出たのはある意味当然かもしれません。
でも結局『なんだかんだやっても予測はできない』で終わらせるのは正直微妙なので、ちょっとあるデータを加工して再び解析を行ってみたら、今度はけっこうおもしろい結果が出るようになりました。でももうこれ以上書いてしまうと長くなってくるので、続きは次回。

解析しているUbuntu(8.04 Hardy Haron)のスクリーンショット。使いやすくて綺麗。

*1 ただ本当にそうかと言われれば正直微妙なところです。株価は人間の心理によって動くわけですから、当然非合理的な動きをする部分もあると思います(たとえそれがマクロ的には消えてしまうにしても)。

Lakai - Fully Flared Intro

Lakai- Fully Flared Intro

Motiongrapherを見ていたらおもしろい動画を発見したので掲載します。最初はただのスケートボードプレーの動画かなぁと思ったら、1分あたりで度肝を抜かれました。なんていうか、無茶苦茶です。

iPhone meets Softbank mobile

この度、ソフトバンクモバイル株式会社は、今年中に日本国内において「iPhone」を発売することにつきまして、アップル社と契約を締結したことを発表いたします。
「iPhone」について | ソフトバンクモバイル株式会社

おお〜！！！！ついに日本でiPhoneが販売されるのか〜！でもまさかソフトバンクとは思わなかった…

正直パケット代とかCMとかどーなるんでしょー？あのお父さんCMで出すんでしょーか？うーん…

Ubuntu 8.04 (Hardy Heron)にRpyをインストールする。

pythonの上からRを扱いたかったので、RPyをUbuntu 8.04(Hardy Heron)にインストールしようと試みたんですが、どうやらこのまま
$ sudo apt-get install python-rpy
とかでインストールすると、インポートの段階でRPyモジュールが足りないなどのエラーがでるようです。

えーそれじゃあソースからコンパイルするしかないのかなぁ面倒くさいなぁとか思っていろいろ調べてみたら、どうやら8.04だけではなく、Ubuntu全般で起こる問題なようでいろんな人が書き込みを行っていました→(rpy import error at Gusty)。自分の他に躓いた人もいるかもしれないのでメモしておきます。

いろんな方法がありますが、いまのとこ一番楽な解決方法は以下のとおり。

1. "http://cran.r-project.org/bin/linux/ubuntu/"にアクセス。自分が使っているヴァージョンを選択(自分の場合は"hardy")
2. "r-base-core_2.7.0-1hardy0_i386.deb"と"python-rpy_1.0.3-1hardy0_i386.deb"の両方をダウンロードしてインストールする(ただし自分の場合AMD64を使っていたので、amd64をインストールしました)。警告がでるかもしれないけど気にしない。
3. 無事インストールが完了したら、
   $ python
   >>> import rpy
   できちんと動作するか確認する。

本来こんなどーでもいいことで時間が取られるのはあんま好きじゃないんですが、いやぁーとりあえず解決できてよかったよかった。

ランダムフォレスト(Random Forest)法を用いて株価予測を行ってみた。-02

第二回。今回は予測を行うにあたって、統計ソフトであるRの"randomForest"パッケージを利用する方法を用いました。

randomForestパッケージのインストール

randomForestパッケージをRに追加するには、CRANからネットワークインストールを行うやり方が最もてっとり早いです。具体的にはRを起動させてから

> install.packages("randomForest")

と入力してください。コンパイルが開始されます。
もしコンパイラが存在しないとかでエラーがでた場合、apt-getかなんかで適当にインストールしてから、再度試みるとうまくいくはずです。

randomForestパッケージの利用

randomForestパッケージの利用手順は、おおまかに以下のようになります。

1. library(randomForest)でパッケージを読み込む。
2. dat.rf <- randomForest(formula ,data=dat ,ntree=500)
   でRFに学習させる。dataには学習に用いるデータテーブルを指定。formulaには推測に用いる変数を"y~x1+x2"のようにして記述する。y以外のすべての変数を推測に用いるならば"y~."とする。ntreeには発生させる決定木の数を指定する。
3. dat.pred <- predict(dat.rf ,newdata=dat.test)
   でRFに推測させる。newdataには推測に用いるデータテーブルを指定する。

※補足

• plot(dat.rf)とすると、決定木の数と予測精度の推移をグラフで見ることができる。
• varImpPlot(dat.rf)とすると、各説明変数の予測に対する寄与度を見ることができる。

詳細は「randomForestパッケージのリファレンス(pdf)」を参照してください。
次回はこのパッケージを用いて実際に予測してみます。

Pythonでiniファイルとかの設定ファイルを読み込んでみる

Python上でiniファイルなどの設定ファイルを読み込むための方法。Pythonは"Battery Included"の発想があるため、特に何も用意しなくても設定ファイルを読み込むライブラリ"ConfigParser"が標準でついてくる。このConfigParserがとっても便利で感動したのでおもわず投稿。

ConfigParserを用いて設定ファイルを読み込むには、ただ以下の数行を付け加えるだけでよい。

import ConfigParser

config = ConfigParser.RawConfigParser()
config.read('settings.conf')
directory = config.items('Directory') #Directoryセクションの読み込み
print directory # [('item1','value1'),('item2','value2'),...] 
print config.get('Directory','item1') # 'value1'

実質設定ファイルの読み込みが3行ですんでる！すごい！

これをC#で行うとなると結構大変で、StringBuilderなどを使って構文解析をする必要がある。
実際はもうiniファイルではなく、XMLなどを用いて設定ファイルを読み込むのが主流になっているはずだからフェアな比較とは言えないけど、ちょっとしたアプリケーションを作るときにXMLは少し気合が入りすぎというか、そこまでしなくても・・・という感じはある。

なのでこうやってたった3行、しかもリファレンスを見ただけで動作がきっちり理解できるpythonはやっぱすごいなと思ったのでした。

Random Survival Forestなるものがあるらしい

(Photo by Keith Marshall , "Savernake Forest Tree")

• CRAN - Package randomSurvivalForest
• Mr.Hemant Ishwaran

Random Forest(RF)自体最近出来た新しい手法ですが、それを発展させたRandom Survival Forest(RSF)なる手法が最近発表されたらしいです。発表されたのは2007年ぐらいで、2008年の4月8日にRのCRANパッケージに追加されました。

ざっと見た感じだと、ブートストラップを用いてN組のデータセットを生成するところまでは同じですが、その後に不必要な枝を剪定しているのが違いでしょうか。そんな劇的にパフォーマンスが上昇するとまではいかなそうですが、RFよりよい精度で推測が行えるかもしれません。ちょっと試してみたいと思います。

てゆーか発展速度早いって！2008年の4月ってつい最近のことですよ。ほんと機械学習の領域はぐんぐん進んでいきますね。

5/28 追記
ちゃんと読んでみたら、このRandomSurvivalForestはRandomForestを生存分析に適応した手法であることがわかりました。なので今回はあんまり使えなさそうな気がします。というか名前まんまですよね。はぁ。

ランダムフォレスト(Random Forest)法を用いて株価予測を行ってみた。-01

さて、今回は2001年ごろに登場した、比較的新しい手法と言われる集団学習の一派『ランダムフォレスト(Random Forest)』法を用いて株価の予測を行ってみました。

ランダムフォレストって？

ランダムフォレスト(RF; random forest)法は,集団学習法の一種である.集団学習(ensemble learning,アンサンブル学習とも呼ぶ)法は決して精度が高いとは言えない複数の結果を組み合わせ,精度を向上させる方法である.いわば,「三人寄れば文殊の知恵」である.集団学習法の中の代表的な方法としてはバギング法,ブースティング法,ランダムフォレスト法がある.
(p.256 金 明哲,村上 征勝『ランダムフォレスト法による文章の書き手の同定』,2007)

大体のイメージとしては、あんまり精度が良くない予測システムを200〜500個くらい作って組み合わせることで、精度の良い結果を得ようとする感じです。

他の手法に比べて利点は？

そもそも「データを元に何かを推測する」というやり方はとっても応用範囲が広い、しかも「推測」っていうのは人工知能と絡んでくるため、ずっとずっと前から研究が続けられてきて、様々な手法が開発されてきました(e.g. Artificial Neural Network(ANN),Genetic Algorism(GA),Support Vector Machine(SVM),...and more!!)。その中の1つとして集団学習があり、さらにその中の1つとしてランダムフォレストがある感じです。

それじゃあランダムフォレストは他の手法と比べて何か利点はあるの？というと、以下の通りとなります。

• 多くのデータセットを用いることによって、とても正確な分類を行うことができる。
• 非常に多くの説明変数を扱うことができる。
• データマイニングにおける分類問題において、説明変数の重要度を見積もることができる。
  (結構重要。他の手法は主にどの変数によってアウトプットが決定されるのか分からないものがある。)
• 欠損したデータを良い精度で推測できるので、データの大部分が欠損していても正確さを保つことができる。
• 従来の手法に比べて、学習速度が早い。

などなど…(英語版のwikipediaよりいくつかの項目を超訳しました)
ランダムフォレストについて詳しく知りたい方は、『Breiman, Leo (2001). "Random Forests". Machine Learning (pdf)』を参照してください。

とりあえず長くなりそうなので今日は概要まで。次回から具体的な手法に入っていきます。

う〜ん…

最近は映像だけじゃなく、投資についての勉強を両方してる感じなんで、もういい加減このブログのタイトルを変えるか、新しくブログをブランチする必要がでてきました。正直ダサいですしねこのタイトル。むぅ…今はまだ考え中ですけど、数日中に決める予定。

どっちにしろ、自分が書く情報がだれか他の人の役に立てるような記事を書きたいと思っています。映像にしろ金融工学にしろシステムトレードにしろ、そういう情報を日本語で公開しているとこって書籍でもネットでも少ないですし。

まだ自分は大学2年生なのであんま大したことは書けないですが、自分ができる最大限のことを実行していきたいと思います。更新量が他のブログと比べて極端にすくないこのブログですが、1人でも読んでくれる人がいるかぎりまだまだやっていくつもりですので、どうぞよろしくお願いします。

Psyop: Happy Valentine's Day

※結構グロいです！閲覧注意！※

久々にPsyop(NYの映像製作集団)を見ていたら出会いました。会社なのにこういうメッセージ性の強い作品を、しかもなんの目的もなく作るあたりがさすがPsyopだなぁと。

メッセージの内容はよくあるシンプルな主張なのですが、それだけにとても難しい。難しいです。

python + BeautifulSoupでHTML解析を行ってみた。

諸事情により株価について整理されたデータを作ることになったんで、pythonとBeautifulSoupを使ってHTMLスクレイピングを行うことにしました。

最初は「えーpythonってインデントで識別するんでしょー慣れなさそー気持ちわるそー」って思ってたんですが、実際やってみると簡単で数時間後にはきちんと動作するスクリプトが組めたのでびっくり。すいません侮ってました。これは便利！

いやぁーそれにしてもすごいっすよBeautifulSoup。なにしろ名前が綺麗っすよね。混沌としたHTMLのSoupからBeautifulな部分を取り出すぜ！っていう思いが伝わってきて凄くいいです。そしてなによりHTML解析が便利すぎる。他の言語はどーやって行っているのかよく分かりませんが、トップレベルで使いやすいことは確かです。

これがHTMLスクレイピングの主要部分↓

soup = BeautifulSoup(rawdata)      
target = soup("tr",{'bgcolor':'#ffffff'},{'align':'right'})[0]('small') 

こんな気持ち悪い文法で動くのはBeautifulSoupだけ！：)

こんな文初めて書きました。いやぁーすごいっす。

保険商品の値段から入院する確率を算出してみた-02

前回の続きです。それでは計算によって求められた確率を載せたいと思います。

結果

結果は以下のとおりです。保険料の金額は「価格.comの保険カテゴリ」で調べました。また、金融工学では通常ローンの利率を国債などの無リスク資産の年率で計算するので、今回はそれに合わせて日本の国債の年利1.37%を用いて計算しました（左上のボタンで全画面表示にできます）。

求められた確率をグラフで表すと以下のようになります。

これではちょっと分かりにくいので、1年のうちに入院をする確率に補正してみましょう。10年間ずっと健康である確率をQとします。1年ずっと健康である確率Pがそれぞれ独立であり、すべて等しいと仮定すると、

P = Q^(1/10) (^はべき乗を表す)

となるので、これを用いて計算すると、

となりました。40代に入ってから、入院の確率が幾何級数的に増加していますね。

求めた確率は結局何を表しているのか？

前回の記事では2つの保険商品から1つの単純な賭けの問題に変形して確率を求めました。これは商品の値段ともらえる値段の期待値が等しくなる確率、言い換えると、「あなたも会社側も、お互い得も損もしない賭け事となる確率」です。この確率はお互いにとって中立であると言う意味で、金融工学では「リスク中立確率(Risk-neutral Probability)」といいます。

実際はそれだと会社側にとって利益がでないので、本当の確率は会社にとって有利な確率、つまりリスク中立確率よりも低い確率であると推測されます。

問題点

今回の計算では10年健康でいたら1回だけボーナスがもらえるという前提で計算をしたんですが、実際は保険料はそのままで、10年ごとに健康だったらボーナスをもらえるという商品なんです。そのため今回の計算だと、10年以降も健康でいる確率は考慮されていないわけです。これはわざわざ10年後に自分にとってさらに有利になる(保険料は同じなのに、歳を取るごとに入院する確率は増えるため)商品をわざわざ解約していることになるため、今回算出した確率は元からリスク中立でない、高めの確率を算出しているということになります。

で、その10年後の確率も考慮して解くってなるとちょい行列計算とか入って専門的になっちゃうんで、今回はやりません。ひまがあったらやってみます。

保険商品の値段から入院する確率を算出してみた-01

photo by wonderferret (creative commons license)

ちょっくら好奇心が沸いたので、金融工学を駆使して医療保険から入院をする確立を算出してみました。専門家にしか分からないような文章にしたらそもそもこれを書く意味がなくなってしまうので、あんま興味ない人や数学にあんま強くない人でもそれなりに楽しめるような内容を目指します。

それでは保険に関する簡単な説明をしたあと、確率のおーまかな算出方法を示して、それによる結果を提示していきます。結果は次回に載せるので、算出方法は興味ないけど結果はちょっと興味あるかもって人は今回は読まないほうがいいかも。

保険商品って？

保険にはダイヤモンドの輝きもなければ、
パソコンの便利さもありません。
けれど目に見えぬこの商品には、
人間の血が通ってます。

ていうCMが流れたことがありますが、結局のところ保険と言うのは金融商品の一種で、誤解を恐れずにいうならば、「自分が不幸になるほうにお金を賭けて、当たったらお金がもらえる商品」です。これによって不幸になったときの損失を減らし、万が一のときに備えることができるわけです。

このように金融商品を用いることで様々なリスクを回避したり低減しようとする手法を、一般に「リスクヘッジ(risk hedge)」といいます。金融工学はこのリスクヘッジ全般に対する理論体系が確立されているので、今回はこれを用いて入院する確率を算出してみます。

確率を算出するのに用いた商品

今回は確率の計算に「アリコのいっしょうおまかせ入院保険」のNAプランとNCプランを用いました。この保険は月払いで保険料を払い、もし

• 入院したら1日に1万円
• 手術になったらさらに10万円
• もし10年間入院しなかったらボーナスとしてNAプランなら10万円、NCプランなら20万円

がもらえる保険です。アリコ側は入院率をどれくらいで見積もっているんでしょうか？

新しく商品を作ってみる

その答えを出すためには、まずこの確率的に複雑な2つの商品をどうにかして単純な賭けの問題に変形する必要があります。そこでためしに次の例を考えてみます。

あなたはこの「自分が不幸になると得をする商品」を自由に売買できたとします。そこであなたはNAプランを1単位売って、NCプランを1単位買うことにしました。そうするとどういうことが起こるでしょうか？

• あなたはNCプランの保険料からNAプランの保険料を差し引いた分の金額を支払う必要があります。
• もしあなたが10年のうちに1回でも入院や手術になった場合、アリコ側からはNCプランの給付金がもらえますが、同時にあなたはNAプランの給付金を売った相手に支払う義務を負うことになります。つまりあなたの利益は±0ということになります。
• もしあなたが10年の間ずっと健康でいるならば、あなたは20万-10万=10万の利益を得ることになります。

これはつまり、2つの商品から「あなたが10年間健康でいれば10万の利益が得られる商品」を作ったことになります。この2つの商品の組み合わせで出来た商品を便宜上『健康プラン』と名づけることにします。

そしてこの健康プランは毎月同じ額の料金を支払う必要があります。これは毎月同額の料金を支払うという点で、車や住宅のローンと同じです。つまりこれは言い換えると「あなたはこの健康プランという商品を分割払いで買っている」ことになります。

いよいよ確率の算出

ここまでくればもうちょいです。あとはこの分割払いで買っている健康プランの値段はいったいどれくらいなのかを算出すればいいわけです。仮にこの健康プランの値段が8万円だったとします。健康だったときにもらえる金額は10万円ですので、この8万円を「10年間ずっと健康である確率の期待値」と見ると、その確率は

10年間ずっと健康である確率 = 8万円 ÷ 10万円 = 80%

となるため、逆に10年間のうち1回でも入院する確率は

100% - 80% = 20%

となるわけです。

長くなったので次回はその計算結果を載せたいと思います。ここまで読んでくれた人ありがとー！

投資のプロははたしてサルに負けるのか？

ちょっと挑戦的なタイトル。現在「金融工学入門」の息抜きとして、投資するまえにこれを見ろ！とクオンツから評判の「ウォール街のランダム・ウォーク」を読みすすめています。まだ読んでいる最中なんですけど、内容的には「バブルに乗って急上昇の株を買う」とか「数日で運用金額を数倍にしようとして、値動きが激しい株を買う」なんていうことはせずに、過去の教訓を学び、金融工学の立場からランダムウォーク理論やファンダメンタル分析を頭に入れて、リスクを考慮した投資をしましょうみたいな感じです。

ちなみに「ランダムウォーク理論」っていうのは端的に言うと、

「微小時間の値動きがブラウン運動となる理論」

もっと分かりやすくいうと、

「すべての値動きはランダムだから過去の値動きから将来の株価の分析やっても意味ないよ！」

っていう、金融工学から出てきた理論です。この理論によると、「目隠しをしたサルに新聞の相場欄めがけてダーツを投げさせ選んだ銘柄でポートフォリオを組んでも、専門家が注意深く選んだポートフォリオとさほど変わりない運用結果が得られる(p.4)」ことになります。

これは正直プロのチャーチストから見れば不快極まりない理論なわけで、もちろん批判も多いのです。個人的には長いスパンでは成り立つとは思いますが、短期ではあんま成り立たない理論だと思います。前の記事であんなこと書いちゃいましたし。

で、話がそれちゃったんですけど、面白いなと思ったのが1960年代あたりにアメリカで起こった「エレクトロニクス・バブル」の紹介です。そのころはとにかく「成長」という言葉がをアメリカ全土をとりまき、特にテクノロジー関係の銘柄はその高い成長率から、非常に高い株価がつけられました。そのため、テクノロジーとは全く関係のない会社も「エレクトロニクス」という名前をつければ株価が急伸するという異常な事態に見舞われました。

社名こそが鍵であった。「トロン」と名のつく会社は、アストロン、ドゥトロン、バルカトロン、トランジトロンなどなど多数であった。また、「ニックス」のつく会社も、サーキトロニクス、スプロニクス、ビデオトロニクスなどがあり、さらに何社か○○エレクトロソニクスという会社もあった。そのきわめつけは、全ての要素を盛り込んだ「パワートロン･ウルトラソニックス」という会社であった。

投資信託業界のリーダーの一つのドレイファス・アンド・カンパニーのジャック・ドレイファスは、この狂気の流行について次のようにコメントした。

「これまで40年間、靴の紐を製造してきた堅実な小規模会社で、株価収益率6倍の会社があったとしよう。この会社の名前を"Shoelace inc."から"Electronics and Silicon Furth-Burners"に変えたとしよう。今日の市場では"Electronics"と"Silicon"の組み合わせは株価収益率15倍に値する。しかし、本当の鍵は"Furth-Burners"という、誰も理解できない言葉に隠されている。誰も分からない言葉と言うのは、総合評価を倍増させる効果がある。ということは、靴紐製造の事業そのものの株価収益率が6倍、"Electronics"と"Silicon"の名前で15倍、合わせて21倍の価値がある。これがさらに"Furth-Burners"という名前のために倍になるため、この新社名の株価収益率は42倍になるのである。」(p.40-41)

完全にバカにしています。ただ、この馬鹿げたバブルを教訓にしたかというとそうでもなく、やっぱり現在でもテクノロジー系統の株価というのはニュースに敏感に反応する傾向があるようです(ライブドアショックとかもありましたし)。テクノロジー関係について良く分かっていない投資家が多いため、そのぶん敏感にならざるを得ないのかもしれません。

もし狙うとしたらニュースがないためにしぼんで割安になっている株を購入して、プラスのニュースが出たときに売るっていう感じでしょうか。あるいは狂気の沙汰となっている状態に空売りして、急落したときに買い戻すとか(諸刃の剣)。

システムトレードの研究-2

第一回の続きです。今回は初めてのプログラムと言うことなので、平均-分散(mean-variance)理論を利用して簡単なプログラムを組んでみます。

プログラムに用いる言語は何が最適？

1. MetaTraderなどの独自言語
   単純なプログラムの場合はこれが一番効率的だが、いちいち覚えなきゃいけないのがめんどくさい。あと今回扱うデータは30万×6通貨と大量なんで、速度や安定性も心配。
2. C or C++
   速度から見れば一番最適なんだけど、コードの行数が単純に多くなる。ポインタ関係のバグとかで手間取るのは勘弁。
3. Rとかの統計言語
   データ解析で非常に活躍するが、バックテストを行うのには不向き？(本来の分野とはかけ離れている)
4. C# , VB.net などの .net 言語
   一番最適なような気がする。システム部分をいちいち製作しなきゃいけないけれど、その分柔軟性があり、なによりC# 3.0からLinqとラムダ式という超強力な武器が手に入る。速度の問題も気になるところはunsafe or 外部dllという手段があるし。
5. Haskell or Ocaml
   関数型言語で実装したほうがコード数が少なくなる場合が多々あるし、なにしろこのタイプの言語は数学系統に強い。現在は見送るが、将来的に学びたい言語。

といろいろ検討した結果、C#となりました。

始めに作ったプログラム

始めに作った売買プログラムは非常に簡単で、以下の手順だけでシステムを組みました。

なお今回取引する金融商品はFX,売買する通貨はオーストラリアドル、スイスフラン、ユーロ、英ポンド、ニュージーランドドル、アメリカドル / 円の6つで、実際の取引に生じるスプレッド、スワップポイントは今回無視しました。

1. 過去10分前までのデータから収益率r、その標準偏差σを計算し、評価関数
   
   を用いてランク付けする。
2. 一番上のランクになった通貨を売買。期待収益率が正ならば対象通貨を買い、負ならば対象通貨を売る。
3. その10分後に対象銘柄を売却か、購入する。

これを1年間に渡って繰り返した場合、100万円の元出がいくらかになったかを計算します。評価に用いた為替レートの1分足データはForexiteから取得しました。

で、作ったプログラムが以下。正直適当です。

問題発生

ところがこのようなシステムを実際に稼動してみると、どうも収益が逆に-20,30万くらい少なくなっている傾向があるようです。これでは良いシステムとは言えません。

そこでそれなら逆転の発想で、単純に売買ルールを逆にしてみればプラスになるんじゃ！ということでシステムを稼動してみると、思惑通りプラスとなりました。今回のバックテストで得られた結果のグラフがこちらです。

カーブの大小はありますが、どれも上方向を向いており、マイナスの結果になっているグラフはありません(ただし月単位で見た場合、負の収益を出している月もある点に注意)。また収益もスプレッドを無視しているとはいえ、含み益が150万～380万。年利換算で150% ～ 380%。平均の月利は約8.0～11.8%という結果になりました。これは明らかに今回の取引が「劣マルチンゲール(有利な取引)」であることを表しています。

なぜそんなことが起こるのか

もし微小時間の値動きがブラウン運動のような完全にランダムな値動きをするならば、このように4年連続で、しかも比較的綺麗な損益曲線を描くことにはならないはずです。このことから以下のように結論付けることにしました。

「少なくとも過去4年間、外国為替市場における為替のレートは微小時間で正または負の方向に値が進んでいる場合、後にその動きとは逆の値動きになる傾向がある。」

なぜこのような値動きになるのか。これはあくまで仮定ですが、自分は以下のように考えました。

• 外国為替市場では売りの圧力と買いの圧力がせめぎあって、その理想均衡点が価値を決定する。
• その理想均衡点の動きは緩やかに推移していくため、微小時間では無視することができる。
• もし価値が理想均衡点から明らかに乖離した値をとる場合、市場は価値が理想均衡点となるように逆方向の圧力がかかり、補正される。

これならば上の結論の説明がつきそうです。そして2007年の収益が最も少なかった理由として、

「2007年はサブプライム問題などにより売り買いの圧力が大きく変動したため、理想均衡点自体が大幅に変動し、予測がしにくかったのではないか」

と推論しました。

次回はこの仮定を検証してみますが、それよりもまずは金融工学をもっと勉強しなきゃ･･･（まだ300ページしか読み終わってない）自分が知らないことはまだまだたくさんある。

追記(03/23)

年利換算が間違っていたので修正しました。資産額じゃなくて純利益だったのを忘れてました。だめだこりゃ。