Search on the blog

2015年3月29日日曜日

Lucasの定理

 Lucasの定理というものを知った。大きな数の重複なし組み合わせの値を素数を法として求めるときに使える。

非負整数m, n, および素数pについて以下が成り立つ。

 

ただし、mi, niはそれぞれ以下を満たす。
証明はここが分かりやすい。

この定理を見て思ったのは、ほとんどゼロになるのでは?ということだった。
例えば、p=2の場合を考える。mi < nとなるようなiが存在すると右辺が0になってしまう。ということは、すべてのiについてmi ≥ niである場合のみ右辺は非零となる。この確率は10桁の二進数だと、pow(3/4, 10) ~ 0.0563程度である。

本当か?ということで試してみた。

#include <iostream>

using namespace std;

int comb[1024][1024];

int main(int argc, char **argv) {
    
    comb[0][0] = 1;
    for (int i = 1; i < 1024; i++) {
        comb[i][0] = comb[i][i] = 1;
        for (int j = 1; j < i; j++) {
            comb[i][j] = (comb[i-1][j-1] + comb[i-1][j]) % 2;
        }
    }

    int total = 0;
    int nonzero = 0;

    for (int i = 0; i < 1024; i++) {
        for (int j = 0; j < 1024; j++) {
            ++total;
            nonzero += comb[i][j];
        }
    }

    cout << 1.*nonzero/total << endl;

    return 0;
}
上のプログラムを実行すると、確かに0.0563と出る。

同様に計算すると、一般にp進数d桁に収まる非負整数m, nでcomb(m, n) % pを計算すると、非零になる割合はpow(1/2 + 1/(2p), d)程度だと考えられる。

Lec 10 | MIT 18.02 Multivariable Calculus, Fall 2007

 多変数微積分学の10回目の講義。second derivative testを使った極値判定について。

Copyright Information
Prof. Denis Auroux, MIT 18.02 Multivariable Calculus, Fall 2007
View the complete course at: http://ocw.mit.edu/18-02F07

講義概要
  1. Second Derivative Test
    1. 特殊例: ax2 + bxy + cy2
    2. 一般例
    3. テイラー展開を使った一般例の説明
ヘッセ行列が出てくると思ったが出てこなかった。二変数の場合の説明なので、それで事足りている。1.の特殊例の話は単なる導入部の例題かと思っていたが、3.の説明に繋がっていた。停留点付近の関数値を二次のテイラー展開近似で表すと、一般の関数でも1.の場合に帰着する。

補足
 ヘッセ行列を使った判定法を忘れていたので、調べてみた。(以下講義の内容ではありません)

  • 停留点xにおいてヘッセ行列の固有値はすべて正 ⇒ xは極小
  • 停留点xにおいてヘッセ行列の固有値はすべて負⇒ xは極大
  • 停留点xでヘッセ行列が正負両方の固有値を持つ ⇒ xは鞍点

これは停留点付近の関数を二次のテイラー展開近似すると分かる。二次のテイラー展開はヘッセ行列の二次形式になる。

ヘッセ行列の固有値がすべて正ということは、ヘッセ行列は正定値行列である。よって停留点からどの方向に移動しても関数値は増加する。よってxは極小。

ヘッセ行列の固有値がすべて負の場合は、停留点からどの方向に移動しても関数値は減少する。よってxは極大。

ヘッセ行列の固有値が正負両方ある場合は、ある固有ベクトル方向に移動すると関数値が増加し、別のある固有ベクトル方向に移動すると関数値が減少する。よってxは鞍点。

2015年3月28日土曜日

Lec 9 | MIT 18.02 Multivariable Calculus, Fall 2007

 多変数微積分学の9回目の講義ビデオを見た。


Copyright Information
Prof. Denis Auroux, MIT 18.02 Multivariable Calculus, Fall 2007
View the complete course at: http://ocw.mit.edu/18-02F07

講義概要
  1. 偏微分と接平面
  2. 停留点
    1. 極大
    2. 極小
    3. 鞍点
  3. 最小二乗法
 停留点は英語でstationary pointというのは知っていたが、critical pointとも言うらしい。講義ではcritical pointという呼称を使っている。
 停留点が出てきたらすぐに二次微分に入って、極大/極小の判定をするのかと思いきや平方完成を使って極小の判定をするところまでで終わっている。二次微分を使った判定は次回の講義で、ということらしい。このワンテンポおく感じが何とも絶妙だ。現在の知識だとこういうアプローチが出来るけど、実はもっと汎用的なやり方があるという伏線を張っている。
 最小二乗法では直線での補完を行ったあと、ムーアの法則に話がうつる。ムーアの法則は指数関数で表記できるが、対数を取ることで直線と同様の手法が使えることを紹介している。

Lec 8 | MIT 18.02 Multivariable Calculus, Fall 2007

 第8回目の講義。相変わらず丁寧である。この人の講義が特別丁寧なのか、アメリカの大学の講義はどれも同じくらい丁寧なのか気になるところ。

Copyright Information
Prof. Denis Auroux, MIT 18.02 Multivariable Calculus, Fall 2007
View the complete course at: http://ocw.mit.edu/18-02F07

講義概要
  1. 多変数関数
    1. グラフ化
    2. 等高線
    3. 偏微分
gnuplotで講義で出てきたグラフを描いてみた。ただグラフを書くだけだとつまらないので、視点を変えながらグラフを眺める様子をgifにしてみた。




Lec 7 | MIT 18.02 Multivariable Calculus, Fall 2007

 多変数微積分学の第7回目の講義。翌日テストがあるらしく、これまで学んだことの復習および演習問題の解説が行われる。演習問題はこれまでの講義を聞いていれば理解できるレベル。

Copyright Information
Prof. Denis Auroux, MIT 18.02 Multivariable Calculus, Fall 2007
View the complete course at: http://ocw.mit.edu/18-02F07

 前回までのビデオはほぼ一語一句完全に聞き取れていたが、今回のビデオは一部よく聞き取れない部分があった。講師は”講義をわかりやすく”ということでゆっくり喋ってくれていたのだなと実感するとともに、早口で喋られるとまだまだ自分のリスニング力では聞き取れないのだなと痛感した。

2015年3月26日木曜日

Lec 6 | MIT 18.02 Multivariable Calculus, Fall 2007

 多変数微積分学の第6回目の講義を聞いた。大学の講義らしくなってきた。

Copyright Information
Prof. Denis Auroux, MIT 18.02 Multivariable Calculus, Fall 2007
View the complete course at: http://ocw.mit.edu/18-02F07

講義概要
  1. ベクトルの微分
    1. 速度
    2. 加速度
  2. ベクトルの積分
    1. 弧長
  3. 応用
 ケプラーは、観測から以下の2つの事実をつきとめた。
  • 太陽のまわりを回る惑星の軌道は同一平面上にある
  • 太陽と惑星を結ぶ直線が一定時間内に通過する面積は一定である
 惑星が太陽のまわりを回るのは太陽の重力によるものだとニュートンによって明らかにされた。しかし、ケプラーの時代はそれが未知であった。

 3.の応用では、ケプラーの主張から惑星の加速度ベクトルは太陽と惑星を結ぶベクトルと平行である(重力)ことを導出する。

2015年3月23日月曜日

Lec 5 | MIT 18.02 Multivariable Calculus, Fall 2007

 多変数微積分の5回目の講義を聞いた。まだまだ余裕で理解できるレベル。

Copyright Information
Prof. Denis Auroux, MIT 18.02 Multivariable Calculus, Fall 2007
View the complete course at: http://ocw.mit.edu/18-02F07

講義概要
  1. 直線の媒介変数表示
  2. 直線と平面の交点
  3. サイクロイド
2.の直線と平面の交点は大学入試レベルだと思われるが、面白い問題。
3.のサイクロイドでは、y~0のときの軌跡をテイラー展開で解析するということをやっていた。

その他
  1. 講義開始のときに「昨日は〇〇を学びました」と言っていた。また、授業完了時には、「よい週末を。次の講義は火曜日です。」と言っていた。つまりこの授業は週に1コマではない。日本の大学だと1つの授業は週に1コマというのが普通だけど、そうではないらしい。どうりでペースがゆったりしているわけだ。
  2. 「もう青色のチョークは使わないよ」で笑いが起きたり、サイクロイドのデモで拍手が起きたり、と学生たちの反応がいい。

2015年3月22日日曜日

Lec 4 | MIT 18.02 Multivariable Calculus, Fall 2007

 多変数微分積分学の第4回目の講義を聞いた。


Copyright Information
Prof. Denis Auroux, MIT 18.02 Multivariable Calculus, Fall 2007
View the complete course at: http://ocw.mit.edu/18-02F07

講義の概要
  1. 平面の式
  2. 3変数線形連立方程式の幾何学的な意味
    1. 3平面が点で交わる
    2. 3平面が線で交わる
    3. 3平面が面で交わる
    4. 3平面が交点を持たない
  3. 斉次線形連立方程式Ax = 0の解
    1. 自明な解[0, 0, 0]のみ (det(A) != 0)
    2. 無限の解 (det(A) = 0で3平面の法線ベクトルが同一平面にある)
  4. 斉次線形連立方程式Ax = bの解
    1. 解が1つ定まる(det(A) != 0)
    2. 解が存在しない、または、無数に存在する(det(A) = 0)
その他気づいたこと
  1. 理科系の授業なのに女性が多い。積極的に質問をするのも女性が多い気がする。
  2. 講師がたまにクイズを出す。学生たちは答えるときに番号が書かれた紙を頭の上にあげる。あの紙はアメリカの大学では必須なアイテムなのだろうか。紙を持ってない人もいて、手で数字を示す学生もいる。

Lec 3 | MIT 18.02 Multivariable Calculus, Fall 2007

 多変数微積分学の第3回目の講義を見た。
外積の応用と行列の基礎の話だった。だいたい知ってる内容だった。

3次元空間に3つの点P1, P2, P3がある。これらの点は平面上にある。
もうひとつの点Pが与えられる。Pは上の3点と同一平面上にあるか?

という問題が面白かった。
  • determinantを使ってparallelepipedの体積が0になるかどうかで判定する方法
  • cross productを使ってnormal vectorを出した後、それが平面上のvectorと直交するかどうかで判定する方法
の2つが紹介されていた。

Copyright Information
Prof. Denis Auroux, MIT 18.02 Multivariable Calculus, Fall 2007
View the complete course at: http://ocw.mit.edu/18-02F07

また、講義の終盤で逆行列を求める方法(小さいサイズの行列に対してのみ有効)が紹介されていた。プログラムを作って動かしてみた。4*4の行列まで動作確認済。
#include <cstdio>
#include <vector>

using namespace std;

typedef vector<double> vec;
typedef vector<vec> mat;

/*
 * remove row r and colum c from a matrix
 */
mat removeRC(const mat &A, int r, int c) {
    int n = A.size();

    mat ret;

    for (int i = 0; i < n; i++) {
        if (i == r) 
            continue;

        vec row;
        for (int j = 0; j < n; j++) {
            if (j == c)
                continue;
            row.push_back(A[i][j]);
        }

        ret.push_back(row);
    }

    return ret;
}

/*
 * calculate determinant of a matrix.
 */
double det(const mat &A) {

    int n = A.size();

    if (n == 1)
        return A[0][0];

    double ret = 0;

    for (int i = 0; i < n; i++) {
        double tmp = A[0][i] * det(removeRC(A, 0, i));
        if (i % 2 == 0)
            ret += tmp;
        else
            ret -= tmp;
    }

    return ret;
}

/*
 * make adjoint of a matrix.
 */
mat adj(const mat &A) {
    
    int n = A.size();
    mat ret(n, vec(n));

    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            ret[j][i] = det(removeRC(A, i, j));
            if ((i + j) % 2)
                ret[j][i] = -ret[j][i];
        }
    }

    return ret;
}

/*
 * make an inverse of a matrix.
 */
mat inv(const mat &A) {

    int n = A.size();
    double d = det(A);
    mat ret = adj(A);
    
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            ret[i][j] /= d;

    return ret;
}

/*
 * test use only
 */
void disp(const mat &A) {
    int n = A.size();
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++)
            printf("%.8lf ", A[i][j]);
        puts("");
    }
    puts("");
}

/*
 * test use only
 */
mat multiply(const mat &A, const mat &B) {
    int n = A.size();
    mat ret(n, vec(n));
    
    for (int i = 0; i < n; i++)
        for (int k = 0; k < n; k++)
            for (int j = 0; j < n; j++)
                ret[i][j] += A[i][k] * B[k][j];
    
    return ret;
}

int main(int argc, char **argv) {

    // 2*2 matrix    
    mat A1 = {
        {2, 1},
        {5, 3}
    };
    disp(multiply(A1, inv(A1)));

    // 3*3 matrix
    mat A2 = {
        {1, 2, 5},
        {1, -1, 1},
        {0, 1, 2}
    };
    disp(multiply(A2, inv(A2)));

    // 4*4 matrix
    mat A3 = {
        {0, -1, -1, 0},
        {3, 0, 1, 1},
        {-3, 0, 1, -1},
        {0, 1, 2, 3},
    };
    disp(multiply(A3, inv(A3)));
    
    return 0;
}

2015年3月21日土曜日

Lec 2 | MIT 18.02 Multivariable Calculus, Fall 2007

 Lecture 2を見た。determinant, cross product, area of parallelogram, volume of parallelepipedなど。まだまだ易しい。

Copyright Information
Prof. Denis Auroux, MIT 18.02 Multivariable Calculus, Fall 2007
View the complete course at: http://ocw.mit.edu/18-02F07

この英語の勉強法かなりいい。英語の勉強も出来、大学の数学の復習もできるなんて、なんという効率の良さ。
やはり、「英語自体を勉強するのではなく、英語で好きなことを勉強する。」という姿勢が大事。

その他気づいたことなど。
  • 9面の黒板を使っている。学生がノートを取りやすいように配慮しているのだろう。
  • 講師は適切なタイミングで「質問はないか?」と問いかける。
  • 学生たちの質問が活発。途中で講師は質問攻めに合う。
  • 学生たちは意外と単純な質問をする。MITだからといってみんな天才というわけではないらしい。
  • 単純な質問に対して、講師は親切すぎるくらい丁寧に答える。そんなことも分からないのか?というような態度を決して取らない。
  • lectureの他にrecitationというのがあるらしい。どんな感じか見てみたい。

2015年3月20日金曜日

Lec 1 | MIT 18.02 Multivariable Calculus, Fall 2007

 昼飯食ってるときに見た動画。
MITの多変数微積分学の講義。1回目の講義ということもあって易しい。

英語の勉強も兼ねて、一通り見てみようか。

Copyright Information
Prof. Denis Auroux, MIT 18.02 Multivariable Calculus, Fall 2007
View the complete course at: http://ocw.mit.edu/18-02F07

ドロップダウンリストのオプションをAjaxで変更する

 ドロップダウンリストAで選択した値に応じてドロップダウンリストBのオプションを変更したいという場合。
 例えば、都道府県DDLで東京都を選んだら、市区町村DDLのオプションが選択した都道府県配下の市区町村になるというような機能など。
 古いシステムだとこういうの毎回ポストしてやってる場合が多いです。これをAjaxで作ってみようというただの暇つぶしです。

出来たもの
デモ
補足ですが、僕はにゃんにゃんよりこじまこちゃんが好きです。

ソース
久しぶりにjQuery使った気がする。
<html>
<head>
  <meta charset="utf-8">
  <title>sample</title>
  <script src="https://code.jquery.com/jquery-1.10.2.js"></script>
  <script>

  $(function() {
    $('#groups').change(function() {
      
      // memberドロップダウンリストのクリア
      $('#members').children().remove();

      // Memberドロップダウンリストに選択されたgroupのメンバーを設定
      var team = $('#groups').val(); 
      if (team != "N/A") {
        $.get("api/members/" + team, function(data) {
          $.each(data, function(val, text) {
            $('#members').append($("<option />").val(val).text(text)); 
          });
        }, "json");
      }

    });
  });

  </script>

</head>
<body>
Team &nbsp;<select id="groups" style="width: 100px;">
<option value="N/A"></option>
<option value="teamA">Team A</option>
<option value="teamK">Team K</option>
<option value="teamB">Team B</option>
</select>  
<br/>
Member &nbsp;<select id="members" style="width: 100px;" />
<body>
</html>

2015年3月13日金曜日

元気が出る動画(4)What is Math About?

数学とは何かを教えてくれる動画。
本当の数学のやり方を教えてくれる動画。


Math is not about numbers, math is not about a calculation, math is not about a logic, math is about "情緒." It's the very act of looking inside of your mind and encountering with your own self, your rich inner universe, your "情緒."

↑ この部分、分野は違うけど、2004 DNC Keynoteを彷彿とさせて感動的です。

2015年3月9日月曜日

formulaで大きめのフォントで数式を書く

 ブログに数式を埋め込む場合、formulaというサイトを使っている。
フォーム内にTexのコマンドを入力すると、数式の画像を作ってくれる。

 パワポに数式を貼り付ける場合も、同じ方法が使えると嬉しいのだが、

e^x = \sum_{n=0}^\infty \frac{x^n}{n!}

のように入力すると、


という画像が出来る。これだとフォントが小さすぎて、拡大してパワポに貼りつけたときに画素が粗くなってしまう。

\mbox{\Huge $\displaystyle
e^x = \sum_{n=0}^\infty \frac{x^n}{n!}
$}
のように入力すると、


のようにフォントサイズが大きい画像が出来る。
これだとパワポに貼りつけても綺麗。

2015年3月8日日曜日

numpy.poly1dと母関数とサイコロ

 「部屋とYシャツと私」みたいなタイトルになってしまいました。

 サイコロを振ったときの確率・統計量を母関数を使って計算するテクニックを知ったのでその紹介をします。
 母関数を扱う際に多項式の計算ライブラリがあると嬉しいです。今回はPythonのnumpy.poly1dクラスを使います。

numpy.poly1dクラスの使用例
多項式の加算、乗算、べき乗、根の計算、微分、特定の点での評価などが出来ます。
from numpy import poly1d as poly

f = poly([1,2,3])    # f(x) = x^2  + 2x + 3
g = poly([2, 1])     # g(x) = 2x + 1

# addition, multiplication and power
print f + g          
print f * g
print g ** 5

# roots of (x^2 - x - 1)
print poly([1, -1, -1]).r

# derivatives
print f.deriv()     # 1st derivative of f 
print f.deriv(2)    # 2nd derivative of f

# value at a value
print f(1)          # f(1) = 1^2 + 2 * 1 + 3
以下が実行結果です。指数の部分とそれ以外の部分を2行に分けて、いい感じに出力してくれます。
kenjih$ python sample.py 
   2
1 x + 4 x + 4
   3     2
2 x + 5 x + 8 x + 3
    5      4      3      2
32 x + 80 x + 80 x + 40 x + 10 x + 1
[ 1.61803399 -0.61803399]
 
2 x + 2
 
2
6

母関数とは?
応用に入る前に母関数について簡単に説明しておきます。

母関数とは、数列{an}に対して
 f(x) = ∑ an xn
と定義される多項式のことです。

f(x)のxnの係数がanと等しいというのがポイントで、うまく利用すると「分割数の計算」や「漸化式の一般項の計算」が出来たりします。

サイコロの統計量と母関数
6つの面を持つサイコロを振ります。
各面が出る確率は等しいとします。

f(x) = (1/6) x + (1/6) x2 + (1/6) x3 + (1/6) x4 + (1/6) x5 + (1/6) x6

という関数を定義します。xnの係数は「サイコロを振った時にnが出る確率」と同じです。

f(x)を1回微分してみます。

f '(x) = (1/6) + (1/6) 2x + (1/6) 3x2 + (1/6) 4x3 + (1/6) 5x4 + (1/6) 6x5

1回微分することで、xnの指数部のnが係数のところに出てきます。これで何が嬉しいかというと、

f '(1) = (1/6) + (1/6) * 2 + (1/6) * 3 + (1/6) * 4 + (1/6) * 5 + (1/6) * 6

となり、x=1を代入すると平均値を計算することが出来ます。

次にf(x)を2回微分してみます。

f ''(x) = 0 + (1/6) * 2 + (1/6) 6x + (1/6) 12x2 + (1/6) 20x3 + (1/6) 30x4

2回微分するxnの係数はn(n-1)になります。1回微分して出てきたnと足すと、n2が出てきます。
つまり、

f '(1) + f ''(1) = (1/6) + (1/6) * 22 + (1/6) * 32 + (1/6) * 42 + (1/6) * 52 + (1/6) * 62

となり、二乗の平均が出せます。よって、

f '(1) + f ''(1) - (f '(1))2

で「二乗の平均 - 平均の二乗」、つまり、分散を計算することが出来ます。

ここまではサイコロを1回ふったときの話を考えました。
次はサイコロを複数回ふった場合、母関数がどうなるかを考えてみます。

f(x)の2乗

((1/6) x + (1/6) x2 + (1/6) x3 + (1/6) x4 + (1/6) x5 + (1/6) x6) 2

を考えてみます。2乗を展開すると、x4 の係数は3/36になることが分かります。
これは、

(1/6) x1 と (1/6) x3の積
(1/6) x2 と (1/6) x2の積
(1/6) x3 と (1/6) x1の積

の和です。

次のように言い換えることも出来ます。

1回目の試行で1が出て、2回目の試行で3が出た確率
1回目の試行で2が出て、2回目の試行で2が出た確率
1回目の試行で3が出て、2回目の試行で1が出た確率

の和。これは「サイコロを2回振った時の目の和が4になる確率」です。

よってf(x)の2乗を計算してxnの係数を見ることで、「サイコロを2回投げたときに目の和がnになる確率」を求めることが出来ます。

さらに一般化すると、f(x)のm乗を計算することで、「サイコロをm回投げた時の目の和の分布」を計算できるということになります。また、f(x)のm乗の微分を行うことで、平均や分散などを計算することができます。

例題
6面のサイコロがある。
サイコロには細工がされていて1,2,3,4,5,6が出る確率が1:2:3:4:5:6となっている。
サイコロを30回振った時に出た目の合計が120以上になる確率を求めよ。

まず、普通に動的計画法で解いてみます。
#include <iostream>

using namespace std;

const int N = 30;  // 30 iterations
const int K = 6;   // 6-sided dice

double dp[N+1][N*K+1];

int sumTo(int k) {
    return k * (k + 1) / 2;
}

int main(int argc, char **argv) {
    
    dp[0][0] = 1.0;
    
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N*K; j++) {
            if (dp[i][j] == 0)
                continue;
            for (int k = 1; k <= K; k++) {
                dp[i+1][j+k] += dp[i][j] * k / sumTo(K);
            }
        }
    }

    double ret = 0;
    for (int i = 120; i <= N*K; i++)
        ret += dp[N][i];

    cout << ret << endl;

    return 0;
}

答えは、0.898544となりました。

次に、母関数を用いた解法。
from numpy import poly1d as poly

dice = range(7)
total = sum(dice)
prob = [1. * x / total for x in dice]

f = poly(prob[::-1])
f = f ** 30

p = f.c[::-1]
print sum(p[120:])

短っっ。実行結果は動的計画で計算したものと一致していました。

ついでに母関数の微分を使って平均、分散を出してみます。
from numpy import poly1d as poly

dice = range(7)
total = sum(dice)
prob = [1. * x / total for x in dice]

f = poly(prob[::-1])
f = f ** 30

mean = f.deriv()(1)
sigma = f.deriv()(1) + f.deriv(2)(1) - mean ** 2

print mean
print sigma

因みに今回の問題の場合は各試行が独立なので、平均・分散は以下のように普通に計算できます。結果は母関数から導出したものと同じになりました。
dice = range(1, 7)
tot = sum(dice)

prob = [1. * x / tot for x in dice]

mean = sum([(x + 1) * p for x, p in enumerate(prob)])
sigma = sum([(x + 1) * (x + 1) * p for x, p in enumerate(prob)]) - mean**2

print mean * 30
print sigma * 30

2015年3月5日木曜日

エンジニア日記(4)こんなとこで働きたい

最近、職場の居心地が良い。

居心地は良いのだけれど、楽しくはない。
ココロオドラないのである。

多分こんなところで働くとココロオドルのではないだろうか、と思う職場像を考えてみた。

1.プログラム好きな人がたくさんいる。
2.アジャイルな開発手法を取り入れている。
3.ドキュメントよりもソースが重視される。
4.テストが自動化されている。
5.ROWEを取り入れている。


2015年3月1日日曜日

おいら、レベル6に上がったっぺ

こんばんはごわす。
おいら、レベル6にあがったっぺ。