2018年2月14日水曜日

Why We Created Julia 何故我々はJuliaを作ったか

Why We Created Julia 何故我々はJuliaを作ったか
14 February, 2012

In short, because we are greedy.
端的に言えば,我々は欲張りだからだ.

We are power Matlab users. Some of us are Lisp hackers. Some are Pythonistas, others Rubyists, still others Perl hackers. There are those of us who used Mathematica before we could grow facial hair. There are those who still can’t grow facial hair. We’ve generated more R plots than any sane person should. C is our desert island programming language.
我々はMatlabのパワーユーザーだ.LispハッカーやPython使いやRuby使いもPythonハッカーもいる.髭が生える前からMathematicaを使っていたのもいるし,未だに髭が生えてない仲間もいる.常識的な人にはオススメしないくらい多くのグラフをR言語で描いてきた。そしてC言語は僕らのユートピアだ.

We love all of these languages; they are wonderful and powerful. For the work we do — scientific computing, machine learning, data mining, large-scale linear algebra, distributed and parallel computing — each one is perfect for some aspects of the work and terrible for others. Each one is a trade-off.
いま挙げた言語は大好きだ,どれも素晴らしいしパワフルだけど,科学計算,機械学習、データ・マイニング,大規模な線形代数演算,分散・並行コンピューティング,といった僕らがやるようなものにはどれも一長一短で,仕事に完璧にはまる機能もあれば何とも使い物にならないものもある.どれもトレード・オフなんだ.
We are greedy: we want more.
僕らは欲張りだ.これじゃ十分じゃない.
We want a language that’s open source, with a liberal license. We want the speed of C with the dynamism of Ruby. We want a language that’s homoiconic, with true macros like Lisp, but with obvious, familiar mathematical notation like Matlab. We want something as usable for general programming as Python, as easy for statistics as R, as natural for string processing as Perl, as powerful for linear algebra as Matlab, as good at gluing programs together as the shell. Something that is dirt simple to learn, yet keeps the most serious hackers happy. We want it interactive and we want it compiled.
僕らが欲しい言語はこんな感じだ.まず、ゆるいライセンスのオープンソースで,Cの速度とRubyの動的さが欲しい.Lispのような真のマクロが使える同図象性のある言語で,Matlabのように分かりやすい数学の記述をしたい.Pythonのように汎用的に使いたいし,Rの統計処理、Perlの文字列処理,Matlabの線形代数計算も要る.シェルのように簡単にいくつかのパーツをつなぎ合わせたい.チョー簡単に習えて、超上級ハッカーも満足する言語.インタラクティブに使えて,かつコンパイルできる言語が欲しい.
(Did we mention it should be as fast as C?)
(そういえば、C言語の実行速度が必要だってのは言ったっけ?)
While we’re being demanding, we want something that provides the distributed power of Hadoop — without the kilobytes of boilerplate Java and XML; without being forced to sift through gigabytes of log files on hundreds of machines to find our bugs. We want the power without the layers of impenetrable complexity. We want to write simple scalar loops that compile down to tight machine code using just the registers on a single CPU. We want to write A*B and launch a thousand computations on a thousand machines, calculating a vast matrix product together.
こんなにもワガママを言った上だけど、Hadoopみたいな大規模分散コンピューティングもやりたい。もちろん、JavaとXMLで何キロバイトも常套句を書きたくないし、数千台のマシンに分散した何ギガバイトものログファイルを読んでデバッグするなんて論外だ。幾層にも重なった複雑さを押しつけられるようなことなく、純粋なパワーが欲しい。単純なスカラーのループを書いたら、一台のCPUのレジスターだけをブン回す機械語のコードが生成されて欲しい。A*Bと書くだけで千の計算をそれぞれ千のマシンに分散して実行して、巨大な行列の積をポンと計算してもらいたい。
We never want to mention types when we don’t feel like it. But when we need polymorphic functions, we want to use generic programming to write an algorithm just once and apply it to an infinite lattice of types; we want to use multiple dispatch to efficiently pick the best method for all of a function’s arguments, from dozens of method definitions, providing common functionality across drastically different types. Despite all this power, we want the language to be simple and clean.
型だって必要ないなら指定したくない。もしポリモーフィックな関数が必要な時には、ジェネリックプログラミングを使ってアルゴリズムを一度だけ書いて、あとは全ての型に使いたい。引数の型とかから自動的にメソッドを選択してくれる多重ディスパッチがあって、共通の機能がまったく違った型にも提供できるようにして欲しい。これだけのパワーがありながらも、言語としてシンプルでクリーンなものがいい。
All this doesn’t seem like too much to ask for, does it?
これって、多くを望みすぎてるとは思わないよね?
Even though we recognize that we are inexcusably greedy, we still want to have it all. About two and a half years ago, we set out to create the language of our greed. It’s not complete, but it’s time for a 1.0 release — the language we’ve created is called Julia. It already delivers on 90% of our ungracious demands, and now it needs the ungracious demands of others to shape it further. So, if you are also a greedy, unreasonable, demanding programmer, we want you to give it a try.
僕らがごまかしようのないほど欲張りなのは分かってるけど、それでもぜんぶ欲しいんだ。二年半ほど前、この欲にまみれた言語を作り始めた。まだ完成してないけど、そろそろ1.0のリリースの時期だ。僕らが作った言語の名前はJulia。すでに僕らの無礼な要求に9割方は応えてくれてるけど、ちゃんとした形になるためには僕ら以外の要求も聞かないといけない。だから、君がもし欲張りで理不尽でわがままなプログラマなら、ちょいとこいつを試してもらいたいんだ。

2018年2月6日火曜日

[Python] モンテカルロ法で円周率を求める

モンテカルロ法を用いて,円周率$\pi$を求めるプログラムのメモです.

例として,直径2cm(当然,半径1vm)の円の面積を求めることを考えます.
直径2cmの円は,1辺が2cmの正方形の中にぴったり入るので,円周率$\pi$を使えば面積は
$2 \times 2 \times \pi = 4\pi$
と求められます.

円周率$\pi$ がわからないときにどうすれば円の面積を求められるかと考えると,ある程度の点をランダムに打って,円の中に入った点と外れた点を数えて,その比率を求めてやれば,大体の円の面積が求められます(別に円でなくても構わない).ランダムということを使って,何かを推定する手法のことをモンテカルロ法と言います.

以下のコードを実行してみます.

import random
import math


def monte_carlo(count):
    ok = 0
    for _ in range(count):
        x = random.random()
        y = random.random()
        if (x * x) + (y * y) <= 1:
            ok += 1
    return 4 * ok / count

print('For    100 {0}' .format(monte_carlo(100)))
print('For   1000 {0}' .format(monte_carlo(1000)))
print('For  10000 {0}' .format(monte_carlo(10000)))
print('For 100000 {0}' .format(monte_carlo(100000)))
(このコードは,人工知能はどのようにして 「名人」を超えたのか?に記載されていたものです)
monte_carloという関数の中でx座標とy座標をランダムに生成させてそれが半径1[cm]の円の中に入った場合の点の数を数えます.countは,試行回数で,何点のランダムの点を生成させるかを表す変数です.

上記のコードを実行してみると,試行回数(ランダムな点を生成させる回数)を増やすにつれて円周率$\pi = 3.141592...$に近づいていくことを確かめることができます.
$ python3 monte_carlo.py
For    100 3.2
For   1000 3.172
For  10000 3.142
For 100000 3.14368

当然,REPLで上記のプログラムを実行しても同じような(というのは乱数 - random - を用いているので全く同じにはならない)結果が得られるはずです.
$ python3
Python 3.6.7 (v3.6.7:6ec5cf24b7, Oct 20 2018, 03:02:14) 
[GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.57)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import random
>>> import math
>>> 
>>> 
>>> def monte_carlo(count):
...     ok = 0
...     for _ in range(count):
...             x = random.random()
...             y = random.random()
...             if (x * x) + (y * y) <= 1:
...                     ok += 1
...     return 4 * ok / count
... 
>>> print('For    100 {0}' .format(monte_carlo(100)))
For    100 3.12
>>> print('For   1000 {0}' .format(monte_carlo(1000)))
For   1000 3.124
>>> print('For  10000 {0}' .format(monte_carlo(10000)))
For  10000 3.1204
>>> print('For 100000 {0}' .format(monte_carlo(100000)))
For 100000 3.14728
>>> 

上記のコードをJupyter Notebookで実行した結果はGitHubで見ることができます.

2018年2月5日月曜日

[Pyrhon] コイントスで表か裏が10回連続で出る確率

エレガントな解法、エレファントな解法 〜モンテカルロ法を添えて〜」という山本一成さんのnoteにあったプログラムを実行してみた際のメモです.

問題は
 コインを100回投げて、表か裏が10回連続で出る確率は?
というものです.
これを実際にコイントスをやってみて確かめるという方法です.
具体的には,コンピュータの中で乱数(ランダム)を発生させ,それでコイントスをしていることとするというものです.

pythonでのプログラムは以下のとおりです.

import numpy

def cointoss(try_count, length):
    x = 0
    prev_ret = None
    for _ in range(try_count):

        ret = numpy.random.randint(0, 2)
        if prev_ret == ret:
            x += 1
            if x >= length:
                return True
        else:
            prev_ret = ret
            x = 1
    return False

def monte_carlo(try_count):
    x = 0.0
    y = 0.0
    for _ in range(try_count):
        if cointoss(100, 10):
            x += 1
        y += 1
    return 100.0 * x / y

print("   100 : {0}%".format(monte_carlo(100)))
print("  1000 : {0}%".format(monte_carlo(1000)))
print(" 10000 : {0}%".format(monte_carlo(10000)))

print("100000 : {0}%".format(monte_carlo(100000)))

これをコンパイル,実行してみると以下のようになりました.
$ python cointoss.py
   100 : 9.0%
  1000 : 9.4%
 10000 : 8.91%
100000 : 8.695%

上記の記事の中では,近似解が
       100回: 14.0%
     1000回 : 8.8%
    10000回: 8.89%
  100000回: 8.773%
となっています.
正解は8.669%だそうです.

ちなみに,プログラムの最後の部分に

print(" 1000000 : {0}%".format(monte_carlo(1000000)))
print("10000000 : {0}%".format(monte_carlo(10000000)))

と書き足して,3回実行させてみた結果は以下のとおりです.

(その1)
$ python cointoss.py
     100 : 10.0%
    1000 : 8.4%
   10000 : 8.94%
  100000 : 8.535%
 1000000 : 8.6777%
10000000 : 8.66218%

(その2)
$ python cointoss.py
     100 : 11.0%
    1000 : 7.5%
   10000 : 8.87%
  100000 : 8.634%
 1000000 : 8.6414%
10000000 : 8.67728%

(その3)
$ python cointoss.py
     100 : 6.0%
    1000 : 10.1%
   10000 : 9.0%
  100000 : 8.815%
 1000000 : 8.5956%
10000000 : 8.68189%