[基]連立一次方程式（第六章）

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


線形代数でやったのをプログラムにしようか！っていうだけ。
だから解説はいらないと思うので関数だけさらしておきます。
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
を十進法表現になおすとどうなるのでしょうか。まず３つの部分にわけてみましょう。
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)
ね？
これで生まれてしまう誤差を丸め誤差と呼びます。


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


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

[基]数値積分と乱数（第六章）

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


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法による積分）

[ス]練習問題（第七章）

授業スライド「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

[基]繰り返しのある関数（第四章）

条件付き繰り返し
三章ですでに繰り返しを扱いましたが、今回も扱います。ただし、条件付き。
そもそも、パソコンのように、根本的には単純なことしかできない機械では「何度も作業を繰り返す」ということは必要不可欠なのです。


さて、ここで「1からnまでの数のうちkの約数の最大を求める」ことを考えましょう。
まあ言い換えれば「nとkの最大公約数を求める」ってことですが、こう考えるとよりややこしくなるのでスルーで。


ではどうプログラムするか考えましょう。
n,n-1,n-2...と順番に、kの約数になるかどうか確かめていけばいいですね。そして、最初に割れた数を答えとすればいい。
でも、for in〜による繰り返しでは、割り切れる数を見つけた時点でストップすることができません。必要なのは、kの約数が見つかるまで、繰り返す関数です。
そこで使うのがwhile文。実際にかいてみましょう。


def  gd_loop(k,n)　 ←関数名をgd_loopとする
  while  k%n != 0      ←「kをnで割ったあまりが0でない」とき
    n = n-1                  ←nを1小さくする
  end                        ←繰り返し命令終了
  n                           ←最終的なn（=kの約数）を出力
end                         ←関数定義終了


while - endは、「kをnで割ったあまりが0でない」限り繰り返されます。すなわち、「kをnで割ったあまりが0」になったら繰り返しが終了し、そのnを出力することになります。




再帰
これも繰り返しの一種です。ただし、少し毛色が違い、また使えるようになるまで少し慣れが必要です。
数学Bの数列でやった漸化式を覚えていますか。あのときは漸化式から数列を求めましたが、今回は、漸化式を作ってそれをそのままプログラムにしてしまおう、というものです。


三章で「1からnの和を求める」関数sum(n)をつくりました。これです。
def  sum(n)
  sum = 0
  for  i  in  1..n
    sum = sum + i
  end
  sum
end


さて、ここで同じ「1からnの和を求める」を漸化式にしてみましょう。簡単ですね。
sum(1)=1 ,  sum(n)=sum(n-1)+n (n=2,3,4...)
です。これがそのまま関数として定義できちゃいます。


def  sum(n)
  if  n>=2                ←「nが2以上」のとき
    sum(n-1) + n     ←「sum(n-1)+n」を出力
  else                      ←「nが2以上でない」すなわち「n=1」のとき
    1                        ←1を出力
  end                       ←条件分岐終了
end

[基]テスト範囲定義まとめ

Rubyを使う前に
☆pwd：今いるフォルダを示す
☆cd ~：ホーム（ターミナルのあるフォルダ）に戻る
☆cd ruby：フォルダ"ruby"に移動する
☆ls：フォルダ内のファイル名を表示する
☆cat bmi.rb：bmi.rbの内容を表示する
 
第一章
☆include(Math)：cosやsqrtなどの数学関数を使う宣言


☆変数名＝式：右辺の式を計算した値を、左辺に書かれた名前の変数にしまう代入命令


☆関数名(式1,…,式n)：式1,…,式nを計算した値を、関数に引き渡す関数呼び出し式


☆def 関数名(変数名1,…,変数名2)
     式
   end
：関数定義。変数名1,…,変数名2について式という関数を定義する。


☆load("ファイル名")：ファイルから関数定義を読み込む命令。




第二章
☆[式0,式1,…,式n]：大きさnの配列を作る。-※


☆式0[式1]：式0が表す配列の式1番目の値を参照する。


☆式0[式1]=式2：式0が表す配列の式1番目の値に式2を代入する。


☆式.length()：式が表す配列の長さを参照する。


☆高次元の数列は、※の式iを数列にする。 


☆show(式)：式が表す配列を画像にして参照する。




第三章
☆if  条件式
     式1
   else
     式2
   end
：条件式が成り立つときは式1を、そうでないときは式2を計算する


☆"文字"：文字列を作る


☆式1+式2：式1、式2が文字列のとき、それらの文字列をつなげた文字列を作る


☆式1[式2..式3]：式1があらわす文字列の式2番目から式3番目までを取り出した文字列を作る


☆for  変数  in  式1..式2
     命令1
      ：
     命令n
   end
：変数の値を式1の値から式2の値まで1ずつ順に変化させながら、命令1から命令nを繰り返し実行する




第四章
☆while  条件式

     命令1

      ：

     命令n

   end
：条件式が成立する間だけ命令1から命令nを繰り返し実行する


☆Array.new(式)：大きさが式の値であるような配列を作る




第五章
配布プログラムbench.rbで定義されている関数の使い方です。


関数の実行時間をグラフにする関数run
☆run("関数名",式)
：関数名(式)の計算を行い、その際の計算時間をX座標が式の位置に表示する


☆run("関数名",式1,式2)
：関数名(式2)の計算を行い、その際の計算時間をX座標が式1の位置に表示する


上記のグラフの表示方法を変更する関数command
☆command("set logscale y")：グラフのY軸を対数スケールにする。X軸も可能


☆command("unset logscale")：対数スケールをやめる


☆command("set xrange [式1:式2]")：X軸の範囲を式1から式2までにする。Y軸も可能


☆command("set autoscale")：表示範囲を自動的に変更する




第六章
☆rand()：0以上1未満の数をランダムで出力する

[ス]練習問題（第五章）

授業スライド「05アルゴリズムと計算量」の練習問題の解答例です。


• 配列 a の(先頭を0 番目としたときの)i 番目 以降の値の中で、最小値が出現する番号を
答える min_index(a,i) を追加して、単純整列 法を完成させよ。


def  min_index(a,i)
  min = i
  for  j  in (i+1)..(a.length()-1)
    if  a[j] < a[min]
      min = j
    end
  end
  min
end




• 省略された部分を補って関数 merge を完成させよ。
（青色が省略されていた部分）


def  merge(a,b)
  c = Array.new(a.length()+b.length())
  ia=0
  ib=0
  ic=0
  while ia < a.length() && ib < b.length()
    if a[ia] < b[ib]
      c[ic]=a[ia]
      ia = ia + 1
    else
      c[ic] = b[ib]
      ib = ib + 1
    end
    ic = ic +1
  end
  if ia < a.length()
    for j in ic..(c.length()-1)
      c[j] = a[ia]
      ia = ia +1
    end
  else
    for k in ic..(c.length()-1)
      c[k] = b[ib]
      ib = ib +1
    end
  end
  c
end

[教]練習問題（第五章）

配布プログラムは、
http://lecture.ecc.u-tokyo.ac.jp/johzu/joho-kagaku/text/code/
からダウンロードしてくださいな。


練習5.1 （再起的アルゴリズムの使用）
略




練習5.2 （正しさの確認）
確認は略。
関数はこちら↓
def  fib_count(k)　　　←補助関数
  j=0
  for  i  in 0 .. k
    if  fibr(i) != fibl(i)
      j=j+1
    end
  end
  j
end


def  fib_check(k) 　　←すべて同じ値ならtrueと返す
  j==0
end




練習5.3 （近似値との比較）
include(Math)
def  fiba(k)
  phi = (1+sqrt(5))/2
  phi**(k+1)/sqrt(5)
end


誤差は自分で調べてね＞＜




練習5.4 （代入の順序）
試してないけど、多分2^(k-1)が出力されるでしょう。
入れ替えたファイルはこうなります。
def  fibl(k)
  f = 1
  p1 = 1
  for  i  in 2..k
    p1 = f
    p2 = p1　　←すなわち、p1=p2=fになるので
    f = p1+p2　←fは、前のfの二倍になる
  end
  f
end
つまり、fは繰り返すごとに二倍になります。




練習5.5 （計算時間の実測）
続く

[基]条件分岐と繰り返し（第三章）

条件分岐
さて、数学をやっていると、よく「場合分け」が必要になりますね。

Rubyで関数をつくるときも、場合分け、すなわち条件分岐が必須になります。

条件式

x>y　　　(>)xがyより大きい

x>=y　　 (≧)xがy以上

x==y　　 (=)xとyが等しい　←「=」じゃなく「==」なのに注意

x<y　　　(<)xがyより小さい

x<=y　　 (≦)xがy以下

x!=y　　  (≠)xとyが異なる

さて、では、xとyのうち大きい値を求める関数を作ってみましょう。

同時に条件分岐の仕方を学びましょう。

def  max2(x,y)　←変数x,yに関する関数、max2(x,y)を定義する

  if  x>=y　　　　←もしx≧yだったら

    x　　　　　 　←xを出力

  else　　　　 　←違う場合（x<yだったら）

    y　　　　　　←yを出力

  end　　　　　 ←条件分岐終了

end　　　　　　←定義終了

条件分岐で使うのは、このif-else-end式です。

さて、もっと複雑な条件分岐をつくる際はどうすればいいのでしょうか。

そんな時は、条件式を組み合わせ、さらにif-else-end式を組み合わせます。

条件式の組み合わせ

x>y || x==0　　x>y または x==0

x<y && y<z　　x<y かつ y<z

!(x<y && y<z)　(x<y かつ y<z)でない

さて、この条件式の組み合わせを使って、

x,y,zのうち一番大きい値を求める関数を作りましょう。

def max3(x,y,z)

  if  x>y && x>z　  ←もし「x>yかつx>z」だったら

    x　　　　　　　  ←xを出力

  else　　　　　　　←違う場合

    if  y>x && y>z　 ←もし「y>xかつy>z」だったら（条件分岐の中に条件分岐）

      y　　　　　　　←yを出力

    else　　　　　　←違う場合（z>xかつz>y）

      z　　　　　　　←zを出力

    end　　　　　　 ←内側の条件分岐終了

  end　　　　　　　←外側の条件分岐終了

end　　　　　　　 ←関数定義終了


できましたか？
少し複雑ですが、練習問題で慣れてください。


論理演算
数学で、「命題」だの「真」だの「偽」だのやった覚えがあると思います。
Rubyでは、変数を含む命題をつくることができます。

たとえば、「xは偶数である」という命題を作りましょう。
作り方は関数を定義するのと同じです。

def  is_even(x)　←命題名
  x%2 == 0　　　←「xが偶数」⇔「x÷2の余りが0」
end

これで、入力したxによって、真すなわちtrue、または偽すなわちfalseを返してくれます。
is_even(2)やis_even(8)とすれば、trueと返しますし、
is_even(5)やis_even(143)とすれば、falseと返します。

これをif文の条件式に使うことも出来ます。
たとえば、「nが偶数だったら1足して、奇数だったらそのまま返す」関数を作るとします。
その場合
def  examle(n)
  if  is_even(n)
    n+1
  else
    n
  end
end
とすればよいのです。


文字列（なんでここにこれwとか思ったら負け）
文字列を値として扱うことができます。「"」で囲むだけです。
s = "abra"
t = "cadabra"
と入力すると、変数sとtの値がそれぞれ"abra"、"cadabra"になります。
そのあと
s+t
と入力するとあら不思議。
"abracadabra"と死の呪文が返ってくるのです。
ちなみに「"」をつければ、数字も文字列として扱えます。

文字列の一部を抜き出すこともできます。
t [1..3]とすると
"aba"と返ってきます。
t[0..0]とすると
"c"と返ってきます。
注意しなければいけないのは、
一文字でも「0から0まで」のように書かなければならないこと。


繰り返し
はいはいはいやってまいりましたーーー
同じ操作を何回も繰り返すってことですよ。
これは練習問題やらないと出来るようにならないです頑張れ。

たとえば、「1からkまでの整数の和」を求める関数sum(k)を考えましょう。
1からkまでを繰り返し足すんですね。


def  sum(k)
  sum = 0　　　　　←和を表す局所関数sum。最初は0。
  for  i  in 1..k　　　←足す数を示す局所変数i。iは1〜kまで繰り返す
    sum = sum + i　 ←iを足していく
  end　　　　　　　←for〜に戻る。繰り返しがkまで終わったら、次の行へ進む。
  sum　　　　　　　←最終的なsumを出力
end


もうお分かりですね…
for (変数)  in  (はじめの数値)..(終わりの数値)
  (命令）
end
が繰り返しの方法です。