『情報科学入門 Rubyを使って学ぶ』の練習問題の解答例の他、授業スライドの練習問題の解答例や、プログラミングのポイント等を掲載します。他クラスの方の利用歓迎。改良・訂正箇所ありましたらご指摘ください。

2011/02/28

[管理情報]アップ状況

試験のために、印刷したいけど、ブログ記事印刷しにくい!って方はご一報ください。
希望のページをPDFにします。


整列法とアライメントはUPしてませんごめんなさい。
スライドみてね。




○→アップ済  △→一部アップ済  ×→未アップ


第一章
 スライド練習問題○・教科書練習問題○・プログラミング基礎○
第二章
 スライド練習問題○・教科書練習問題○・プログラミング基礎○
第三章
 スライド練習問題○・教科書練習問題○・プログラミング基礎○
第四章
 スライド練習問題○・教科書練習問題○・プログラミング基礎○
第五章
 スライド練習問題○・教科書練習問題△(練習5.5まで)・プログラミング基礎△
第六章
 スライド練習問題○・教科書練習問題○・プログラミング基礎○
                              ↑New(2011/01/29)
第七章
 スライド練習問題○・教科書練習問題×・プログラミング基礎×



教科書公式配布プログラムはこちらから
http://lecture.ecc.u-tokyo.ac.jp/johzu/joho-kagaku/text/code/

伊地知先生のスライドはこちらから
http://lecture.ecc.u-tokyo.ac.jp/~cichiji/is-10/is-10.html



質問は、コメントまたは、is.ruby2010☆gmail.com(☆を@に)までメールください。

2011/02/06

[基]アルゴリズムと計算量(第五章)

この章で言いたいことを一言でいうと、
ある問題を特にしても、アルゴリズム(解き方)が違えば、実行時間が違うよ!
ということです。


皆さんも、繰り返しや再帰を使う関数を実行していたとき、変数を大きな数にしすぎてフリーズした!なんて経験があるはず。
あれは実際、コンピュータが超ガンバって計算してるんですよ。
コンピュータだって一瞬で全ての答えを出すわけではありません。


だけど、それはつまり、
同じ答えを導くにしても、計算量が減るようにアルゴリズムを改良してやれば、より素早く答えを出せる
ということ。
教科書では、例として、フィボナッチ数を求めるアルゴリズムと、整列のアルゴリズムを照会しています。


フィボナッチ数のアルゴリズムはこちら。


その1.再帰的アルゴリズムfibr(k)→実行時間はkの指数関数に近似できる
def  firb(k)
  if  k==0 ||  k==1
    1
  else
    fibr(k-1) + fibr(k-2)
  end
end


その2.数え上げアルゴリズムfibl(k)
def  fibl(k)
  f=1
  p1=1
  for  i  in 2..k
    p2=p1
    p1=f
    f=p1+p2
  end
  f
end


その2の2.数え上げでfibの下6桁だけを求めるfibl6(k)
             →実行時間はkの一次関数に近似できる
def  fibl6(k)
  f=1
  p1=1
  for  i  in  2..k
    p2=p1
    p1=f
    f=(p1+p2)%/1000000
  end
  f
end


その3.黄金比による近似解アルゴリズムfiba(k)
def  fiba(k)
  phi=(1+sqrt(5))/2 ←これが黄金比
  phi**(k+1)/sqrt(5)
end





これらのアルゴリズムの速さを比べるために用いられるのが、計算量とそのオーダー。
ですが、これは前期の必修科目情報の範囲なので、とばします。
詳しくは教科書P86〜89をどうぞ。
参考までに、fibr(k)はO(φ^k)、fibl(k)はO(k)です。


さらに!

その4.行列を用いたアルゴリズムfibm(k)


フィボナッチ数を、vk=( fib(k+1) , fib(k) )として考えると、2×2行列Q=[[1,1],[1,0]]をつかって、
 v0=( 1 , 1 )   ,   vk+1=Qvk (k=0,2,3,...) 
すなわち、
  vk=Q^k × v0
と書けます。


行列Qのn乗計算をそのままやるとO(n)になり、数え上げプログラムと大差なくなってしまいますが、
Q = E   (n=0)
      (Q^(n/2))^2  (nは2以上の偶数)
      Q×Q^(n-1)   (nは奇数)
としてやれば、あら不思議、めっちゃ計算量が減っちゃいますね。
なんと、O(logk)になります。


def  matmul(a,b)  ←行列a,bの積を求める関数
  c = make2d(2,2)
  c[0][0] = a[0][0]*b[0][0] + a[0][1]*b[1][0]
  c[0][1] = a[0][0]*b[0][1] + a[0][1]*b[1][1]
  c[1][0] = a[1][0]*b[0][0] + a[1][1]*b[1][0]
  c[1][1] = a[1][0]*b[0][1] + a[1][1]*b[1][1]
  c
end


def  matsquare(a)  ←行列aの自乗を求める関数

  c = make2d(2,2)
  c[0][0] = a[0][0]**2 + a[0][1]*a[1][0]
  c[0][1] = a[0][0]*a[0][1] + a[0][1]*a[1][1]
  c[1][0] = a[1][0]*a[0][0] + a[1][1]*a[1][0]
  c[1][1] = a[1][0]*b[0][1] + a[1][1]**2
  c
end

def  matpower(a,n)  ←行列aのn乗を求める関数
  if  n==0

    [[1,0],[0,1]]

  else

    if  n%2==0

      matsquare(matpower(a,n/2))

    else

      matmul(a,matpower(a,n-1))

    end

  end

end



def  fibm(k)

  q=matpower([[1,1],[1,0]],k-1)

  q[0][0]+q[0][1]

end

2011/01/30

[基]連立一次方程式(第六章)

連立一次方程式を行列で解くプログラム


線形代数でやったのをプログラムにしようか!っていうだけ。
だから解説はいらないと思うので関数だけさらしておきます。
jgp(a):行列aを解く本関数
maxrow(a,k):k列目の係数の絶対値が最大となる行をk行目以降から探す補助関数
abs(x):xの絶対値を出す補助関数
swap(a,i,j):配列aのi項目とj項目を入れ替える補助関数


ただ、誤差が生じるのを防ぐために、k行目使ってk列を消去する前に、|ai,k|が最も大きい行i=pを選び、k行目とp行目を入れ替える必要があるのに注意。
※ak,kが0だったら計算できない。また0に近いと、計算途中に非常に大きな値が大量に発生し、情報落ち誤差が生まれてしまう。


def  gjp(a)
  row = a.length()
  col = a[0].length()
  for  k  in 0..(col-2)
    max = maxrow(a,k)
    swap(a,k,max)
    akk = a[k][k]
    for  i  in 0..(col-1)
      a[k][i] = a[k][i]*1.0 / akk
    end
    for  i  in 0..(row-1)
      if  i != k
        aik = a[i][k]
        for  j  in k..(col-1)
          a[i][j] = a[i][j] - aik*a[k][j]
        end
      end
    end
  end
  a
end


def  maxrow(a,k)
  s = k
  t = abs(a[k][k])
  if  k != a.length()-1
    for  y  in (k+1)..(a.length()-1)
      if  t < abs(a[y][k])
        s = y
        t = abs(a[y][k])
      end
    end
  end
  s
end


def  abs(x)
  if  x<0
    -x
  else
    x
  end
end


def  swap(a,i,j)
  t = a[i]
  a[i] = a[j]
  a[j] =t
  a
end

[基]浮動小数点表現と誤差

ここでは二進法が扱われていますが、二進法と十進法の相互変換はできること前提で記述しています。


浮動小数点数
なんかめんどくさそーですけど、まあ実際面倒です。でも数学より簡単。
コンピュータは、すべてのデータをONとOFF、すなわち1と0で扱っています。二進法ってやつです。
で、大きな数を二進数で表すと、すっごい桁数がいる、っていうのはわかると思います。
そこで、なるべく少ない桁数で大きな数字を表現するために用いるのが、浮動小数点法。


たとえば、数字を表現するのに、32ビット与えられたとします。
普通の二進法で表現すると、最大値はたったの4,294,967,295(=2^32-1)。
しかし、浮動小数点数で表現すれば、なんと最大値は約3.4×10^38!なんというビットの節約。すばらしい。


浮動小数点表現にも何種類かありますが、現在、主に単精度(32ビット列で表す)と倍精度(64ビット列で表す)が使われています。
この二つは、ビット数がちょっと違うだけなので、ここでは単精度で説明していきます。


さて、単精度浮動小数点表現では、32ビットで数字を表しますが、まず、このビット列は3つの部分に分けられています。
それが、符号部s指数部d'1..d'c仮数部d1..dmです。
単精度では、それぞれ1桁、8桁、23桁となっています。
これが
  (-1)^s(2) × 1.d1d2...dm(2) × 2^( d'1d'2..d'c(2) - b )
というふうに数字を表現します。(2)は二進数ということです。
bバイアス値といい単精度では127です。


たとえば、次のような単精度浮動小数点数
1111111101100000000000000000000
を十進法表現になおすとどうなるのでしょうか。まず3つの部分にわけてみましょう。
1  11111110  11000000000000000000000
ということは、十進数では
   (-1)^1(2) × 1.11000000000000000000000(2) × 2^( 11111110(2) - 127 )
= (-1)^1 × { 2^0 + 2^(-1) + 2^(-2) } × 2^( 254 - 127)
= -1.75 × 2^127
≒ -2.977471 × 10^38
となります。


逆に、十進数で表された数、例えば
3145728
を単数度浮動小数点表現にしてみましょう。こちらの方が簡単だと思います。
素因数分解をすればいいだけだからです。
   3145728
= 2^20 × 3
= (-1)^0 × 1.5 × 2^21
= (-1)^0 × { 2^0 + 2^(-1) } × 2^( 148 -127 )
= (-1)^0(2) × 1.10000000000000000000000(2) × 2^( 10010100(2) - 127 )
ですね。というわけで、符号部、指数部、仮数部の順にならべて、
01001010010000000000000000000000
となります。
注意するのは、仮数部(1.d1..dm(2))は、1以上2未満の数字しか表現できないということ。


そして、以上の方法では0を表せないのがわかると思います。
なので、特別に、指数部と仮数部がともに全て0のときは0を表すことになっています。


因みに、倍精度浮動小数点表現では
指数部が11桁、仮数部が52桁、バイアス値は1023です。




誤差
浮動小数点数を用いた計算は、近似値による計算なので、誤差が生じます。
十進数では有限の桁数で表せるのに、二進数では無限に続く少数になってしまうような数が存在するからです。
たとえば0.1もそうです。
0.1 = 2^(-4) + 2^(-5) + 2^(-8) + 2^(-9)......
      = 0.0001100110011......(2)
      = 1.100110011......(2) × 2^(-4)
ね?
これで生まれてしまう誤差を丸め誤差と呼びます。


誤差には他に桁落ち誤差、情報落ち誤差、打ち切り誤差があります。


桁落ち誤差:非常に近い2つの値の差を求めた場合に有効桁数が減ることで引き起こされる誤差
情報落ち誤差:大きな数に小さな数を足した場合に、小さな数が大きな数の有効数字の範囲外になることで無視されることで起きる誤差
打ち切り誤差:無限級数によって表されている値の計算を打ち切ることで生まれる誤差

2011/01/29

[基]数値積分と乱数(第六章)

私は数学が苦手なのです。


Rubyは、我々やMathematicaのような積分のテクニックをもっていません。
なので、台形公式、Simpson公式、Monte Carlo法などで近似的に数値積分をおこないます。
わかっておかなければならないのは、これらは近似にすぎず、誤差が場合によってはかなり生まれてしまうことです。


台形公式















この数値積分を行う関数をtrapezoid(xs,xe,n)とすると、定義はこうなります。
def  trapezoid(xs,xe,n)
  deltax = (xe-xs)*1.0 / n
  sum = (f(xs)+f(xe)) / 2.0
  for  i  in 1..(n-1)
    sum = sum + f(xs+i*deltax)
  end
  deltax * sum
end
f(x)を別に定義すれば、いろんな関数を台形公式で数値積分できるのです。


Simpson公式
こちらは、より正確にするため、直線ではなく二次曲線で近似しようぜ!ってもの。















これをsimpson(xs,xe,n)として定義すると、以下のようになります。
def  simpson(xs,xe,n)
  deltax=(xe-xs)/(2.0*n)
  sum=f(xs)+4*f(xs+deltax)+f(xe)
  for  i  in 1..(n-1)
    sum=sum+2*(xs+2*i*deltax)+4*f(xs+(2*i+1)*deltax)
  end
  deltax*sum/3.0

end


Monte Carlo法
Monte Carlo法は、乱数を用いた数学的問題の解法の総称で、積分以外の問題にも用いられます。
それらの問題は、主に決定論的問題(積分含む)と非決定論的問題にわけられます。
前者は、厳密な解が存在し、乱数でそれを近似的に求めるもの、後者は、自然科学・社会科学におけるシミュレーションのように、確率的な変動を含む、厳密な解が存在しないものを指します。
積分において、Monte Carlo法は、「積分区間、その区間におけるf(x)の値域が、0以上有限値以下におさまっている高次元の関数」に対して使われます。


例として円周率を求めてみましょう。

















まあこうやって、ランダムにn個点をうちまくって、そのうち何割が関数より下(または上)にくるのかな?って感じで、積分する・・・ていうか、面積を出します。
で、上の場合の積分を行う関数montecarlo(n)はこちら。
def  montecarlo(n)
  m=0
  for  i  in 1..n
    x=rand()
    y=rand()
    if  x*x+y*y < 1.0
      m=m+1
    end
  end
  m*1.0/n
end


rand()が乱数で、[0,1)の数をランダムに出力します。
乱数が含まれているので、montecarlo(n)は、nを同じ数にしても、毎回違う値になります。
そして、最終的に円周率を出すにはこうですね
def  pi(n)
  4.0*montecarlo(n)
end


ちなみに、参考までに、使い回しできるmonte carlo法の関数はこちら。
積分区間[xs,xe]、その区間のf(x)の最大値がmaxの場合の関数です。
また、xs≧0かつ、f(x)≧0(xs≦x≦xe)でなくては使えません。
def  montecurlo_super(xs,xe,max,n)
  m = 0
  for  i  in 1..n
    x = rand() * (xe-xs)
    y = rand() * max
    if  y<f(x+xs)
      m=m+1
    end
  end
  m*(xe-xs)*max*1.0/n
end

試験対策オススメ練習問題

授業でやったことを理解してるか確認するのにいいんじゃん?と
個人的に思った練習問題一覧。
まあ各練習問題のa)だけをやるってのもいいかもしれません。


第一章
練習1.2 (関数の定義)
練習1.9 (もっと関数の定義)


第二章
練習2.1 (画像データの作成)
練習2.2 (画像データの操作)
練習2.4 (国旗)


第三章
練習3.2 (繰り返しによる画像の作成)
練習3.3 (条件分岐)
練習3.5 (真偽値を求める関数)
練習3.7 (文字列)
練習3.9 (数列を作る)
練習3.10 (いろいろな図形)※チャレンジ問題


第四章
練習4.1 (素数の判定)
練習4.2 (組み合わせ数)
練習4.6 (繰り返しによって値を集める)
練習4.7 (繰り返しによって条件を満たす値を集める)
練習4.9 (再帰的に値を集める)
練習4.10 (再帰的に条件を満たす値を集める)
練習4.13 (配列の並べ替え)


第六章
練習6.2 (Simpson公式による積分)
練習6.5 (浮動小数点数の表現)
練習6.9 (Pivotingの完成)
練習6.14 (解析的な積分が難しい関数の数値積分)
練習6.15 (Monte Carlo法による積分)

2011/01/28

[ス]練習問題(第七章)

授業スライド「07パターン認識入門」の練習問題の解答例です。


•align.rbを完成させて、ATAGとAACのアラインメントを求めよう。
(青色が省略されていた部分)
def  align(s,t)
  m=s.length()
  n=t.length()
  a=make2d(m+1,n+1)
  for  j  in  1..n
    a[0][j] = a[0][j-1] + g()
  end
  for  i  in  1..m
    a[i][0] = a[i-1][0]+g()
  end
  for  i  in 1..m
    for  j  in  1..n
      x = a[i-1][j] + g()
      y = a[i][j-1] + g()
      z = a[i-1][j-1] + q(s[i-1],t[j-1])
      a[i][j] = max(x,max(y,z))
    end
  end
  a
end





•RNase_P.rbをロードすると、文字列を返す(引数なしの)関数seq0()とseq1()が定義される。これらのアラインメントを求めよう。
やってみそ。241になります(右下の要素)


参考までに、traceback(a,s,t)の完成形
def  traceback(a,s,t)
  u = ""
  v = ""
  i = s.length()
  j = t.length()
  while  i>0 || j>0
    if  j>0 && a[i][j] = a[i][j-1] + g()
      u = "-" + u
      v = t[j-1..j-1] + v
      j = j-1
    else
      if  i>o && j>0 && a[i][j] == a[i-1][j-1] + q(s[i-1],t[j-1])
        u = s[i-1..i-1] + u
        v = t[j-1..j-1] + v
        i = i-1
        j = j-1
      else
        if  i>0 && a[i][j] = a[i-1][j] + g()
          u = s[i-1..i-1] + u
          v = "-" + v
          i = i-1
        end
      end
    end
  end
  [u,v]
end