やろーじだい

ブログです

Aizu SICP 第 13, 14 回

先月末にやりました。練習問題タイムで楽しかった。

練習問題 2.11

区間の両端点の符号をテストすると mul-interval は 9 パターン に場合分けできて、3 回以上のかけ算が必要になるのはその中のひとつだけだよ。”

(Page 101).

参加者により掛ける区間 (区間1) と掛けられる区間 (区間2) の下限と上限が、それぞれ 0 より大きいか小さいかで 9 パターンに分けられることに気付いた。区間1, 区間2, 演算結果の区間をそれぞれ X, Y, Z とし、下限を L, 上限を U で表すと以下のような表で表せる。

0 < X.L < X.U X.L <= 0 <= X.U X.L < X.U < 0
0 < Y.L < Y.U
Y.L <= 0 <= Y.U
Y.L < Y.U < 0

この時のそれぞれの場合の掛け算は、上限と下限について以下の表のように考えることができる。

0 < X.L < X.U X.L <= 0 <= X.U X.L < X.U < 0
0 < Y.L < Y.U Z.L = X.L * Y.L Z.L = X.L * Y.U Z.L = X.L * Y.U
Z.U = X.U * Y.U Z.U = X.U * Y.U Z.U = X.U * Y.L
Y.L <= 0 <= Y.U Z.L = X.U * Y.L Z.L = X.L * Y.U
Z.U = X.U * Y.U Z.U = X.L * Y.L
Y.L < Y.U < 0 Z.L = X.U * Y.L Z.L = X.U * Y.L Z.L = X.U * Y.U
Z.U = X.L * Y.U Z.U = X.L * Y.L Z.U = X.L * Y.L

たとえば両方の区間の下限と上限が全て 0 より大きい場合、結果の区間の下限はそれぞれの区間の下限の乗算の結果になる。結果の区間の上限も同様にそれぞれの区間の上限の乗算の結果になるというのがわかる。負の場合などについても符号に気を付ければ同様に考えることができる。よって表からわかるように3回以上の掛け算が必要となるのは、中央の場合のみである。中央のケースは記号だけではわからないため、以下のような元の mul-interval と同じ計算を行った。

(define (mul-interval x y)
  (let ((p1 (* (lower-bound x) (lower-bound y)))
        (p2 (* (lower-bound x) (upper-bound y)))
        (p3 (* (upper-bound x) (lower-bound y)))
        (p4 (* (upper-bound x) (upper-bound y))))
    (make-interval (min p1 p2 p3 p4)
                   (max p1 p2 p3 p4))))

ただし、xy の下限は負で上限は正であることはわかっているので、下限は正と負の積のどちらかに限られ、上限は符号が一致する積のどちらかに限られる。したがってこの部分の最小限の実装は以下のようになる。

(let ((p1 (* (lower-bound x) (lower-bound y)))
      (p2 (* (lower-bound x) (upper-bound y)))
      (p3 (* (upper-bound x) (lower-bound y)))
      (p4 (* (upper-bound x) (upper-bound y))))
  (make-interval (min p2 p3)
                 (max p1 p4)))

 全体のプログラムについてはとくに綺麗に書く方法を思い付かず、単にそれぞれの区間を場合分けで3種類に分類して場合分けで書いた。

(define (interval-type int)
  (cond [(and (> (lower-bound int) 0) (> (upper-bound int) 0))
         :positive]
        [(and (<= (lower-bound int) 0) (>= (upper-bound int) 0))
         :include-0]
        [(and (< (lower-bound int) 0) (< (upper-bound int) 0))
         :negative]))

(define (new-mul-interval x y)
  (let ([x-type (interval-type x)]
        [y-type (interval-type y)])
    (cond [(equal? x-type :positive)
           (cond [(equal? y-type :positive)
                  (make-interval (* (lower-bound x) (lower-bound y))
                                 (* (upper-bound x) (upper-bound y)))]
                 [(equal? y-type :include-0)
                  (make-interval (* (upper-bound x) (lower-bound y))
                                 (* (upper-bound x) (upper-bound y)))]
                 [(equal? y-type :negative)
                  (make-interval (* (upper-bound x) (lower-bound y))
                                 (* (lower-bound x) (upper-bound y)))])]
          [(equal? x-type :include-0)
           (cond [(equal? y-type :positive)
                  (make-interval (* (lower-bound x) (upper-bound y))
                                 (* (upper-bound x) (upper-bound y)))]
                 [(equal? y-type :include-0)
                  (let ((p1 (* (lower-bound x) (lower-bound y)))
                        (p2 (* (lower-bound x) (upper-bound y)))
                        (p3 (* (upper-bound x) (lower-bound y)))
                        (p4 (* (upper-bound x) (upper-bound y))))
                    (make-interval (min p2 p3)
                                   (max p1 p4)))]
                 [(equal? y-type :negative)
                  (make-interval (* (upper-bound x) (lower-bound y))
                                 (* (lower-bound x) (lower-bound y)))])]
          [(equal? x-type :negative)
           (cond [(equal? y-type :positive)
                  (make-interval (* (lower-bound x) (upper-bound y))
                                 (* (upper-bound x) (lower-bound y)))]
                 [(equal? y-type :include-0)
                  (make-interval (* (lower-bound x) (upper-bound y))
                                 (* (lower-bound x) (lower-bound y)))]
                 [(equal? y-type :negative)
                  (make-interval (* (upper-bound x) (upper-bound y))
                                 (* (lower-bound x) (lower-bound y)))])])))
;; Tests

(test-section "interval-type")

(test* "interval-type1"
       :positive
       (interval-type (make-interval 1 2)))

(test* "interval-type2"
       :include-0
       (interval-type (make-interval -1 1)))

(test* "interval-type3"
       :negative
       (interval-type (make-interval -3 -2)))

(test* "interval-type4"
       :include-0
       (interval-type (make-interval 0 1)))

(test-section "new-mul-interval")
(test-section "x+y+")
(test* "new-mul-interval"
       (mul-interval (make-interval 1 5) (make-interval 1 1))
       (new-mul-interval (make-interval 1 5) (make-interval 1 1)))

(test-section "x+y-")
(test* "new-mul-interval"
       (mul-interval (make-interval 1 5) (make-interval 1 1))
       (new-mul-interval (make-interval 1 5) (make-interval 1 1)))

(test-section "x+y+-")
(test* "new-mul-interval"
       (mul-interval (make-interval 1 5) (make-interval -1 1))
       (new-mul-interval (make-interval 1 5) (make-interval -1 1)))

(test-section "x-y+")
(test* "new-mul-interval"
       (mul-interval (make-interval -5 -1) (make-interval 1 2))
       (new-mul-interval (make-interval -5 -1) (make-interval 1 2)))

(test-section "x-y+-")
(test* "new-mul-interval"
       (mul-interval (make-interval -5 -1) (make-interval -1 2))
       (new-mul-interval (make-interval -5 -1) (make-interval -1 2)))

(test-section "x-y-")
(test* "new-mul-interval"
       (mul-interval (make-interval -5 -1) (make-interval -2 -1))
       (new-mul-interval (make-interval -5 -1) (make-interval -2 -1)))

(test-section "x+-y+")
(test* "new-mul-interval"
       (mul-interval (make-interval -5 5) (make-interval 1 2))
       (new-mul-interval (make-interval -5 5) (make-interval 1 2)))

(test-section "x+-y+-")
(test* "new-mul-interval"
       (mul-interval (make-interval 0 5) (make-interval 0 1))
       (new-mul-interval (make-interval 0 5) (make-interval 0 1)))
(test* "new-mul-interval"
       (mul-interval (make-interval 1 -1) (make-interval 5 -5))
       (new-mul-interval (make-interval 1 -1) (make-interval 5 -5)))

(test-section "x+-y-")
(test* "mul-interval"
       (mul-interval (make-interval -5 5) (make-interval -2 -1))
       (new-mul-interval (make-interval -5 5) (make-interval -2 -1)))

以下ではマクロを使って書いてみたが、単に (positive-interal? interval) などの述語を作ればよいだけだった。interval-type の実装があまり良くなかった。折角なので一応載せておく。

(define-syntax match-interval-type
  (syntax-rules ()
    ((_ target (type1 case1) (type2 case2) (type3 case3))
     (let ([target-type (interval-type target)])
       (cond [(equal? target-type type1)
              case1]
             [(equal? target-type type2)
              case2]
             [(equal? target-type type3)
              case3])))))

(define (new-mul-interval x y)
  (match-interval-type x
    [:positive (match-interval-type y
                 [:positive
                  (make-interval (* (lower-bound x) (lower-bound y))
                                 (* (upper-bound x) (upper-bound y)))]
                 [:include-0
                  (make-interval (* (upper-bound x) (lower-bound y))
                                 (* (upper-bound x) (upper-bound y)))]
                 [:negative
                  (make-interval (* (upper-bound x) (lower-bound y))
                                 (* (lower-bound x) (upper-bound y)))])]
    [:include-0 (match-interval-type y
                  [:positive
                   (make-interval (* (lower-bound x) (upper-bound y))
                                  (* (upper-bound x) (upper-bound y)))]
                  [:include-0
                   (let ((p1 (* (lower-bound x) (lower-bound y)))
                         (p2 (* (lower-bound x) (upper-bound y)))
                         (p3 (* (upper-bound x) (lower-bound y)))
                         (p4 (* (upper-bound x) (upper-bound y))))
                     (make-interval (min p2 p3)
                                    (max p1 p4)))]
                  [:negative
                   (make-interval (* (upper-bound x) (lower-bound y))
                                  (* (lower-bound x) (lower-bound y)))])]
    [:negative (match-interval-type y
                 [:positive
                  (make-interval (* (lower-bound x) (upper-bound y))
                                 (* (upper-bound x) (lower-bound y)))]
                 [:include-0
                  (make-interval (* (lower-bound x) (upper-bound y))
                                 (* (lower-bound x) (lower-bound y)))]
                 [:negative
                  (make-interval (* (upper-bound x) (upper-bound y))
                                 (* (lower-bound x) (lower-bound y)))])]))

本文

(define (width i)
  (/ (- (upper-bound i) (lower-bound i)) 2))

本文の中で定義されている幅 (width) だが、区間の幅というのは直感的には 上限 - 下限 であると思ったのだが (上限 - 下限) / 2 と定義されている。つまり中央値から端点への距離を幅として定義しているので注意したい。

練習問題 2.12

自分は make-interval を使って書いた。しかし、参加者の make-center-width を使ったものの方が綺麗なのでそちらも載せた。

;; make-interval 版
(define (make-center-percent c p)
  (let ([error (/ p 100)])
    (make-interval (- c (* c error)) (+ c (* c error)))))

;; make-center-with 版
(define (make-center-percent c p)
  (make-center-width c (* (* c p) 0.01)))

percent に関しても同様で、自分は何も考えず関数を書いてしまったが、widthcenter を使うべきだった (折角本文で作っているのを無視してしまった)

;; idiot 版
(define (percent i)
  (let* [(c (center i))
         (ub (upper-bound i))]
    (* (- (/ ub c) 1) 100)))

;; genius 版
(define (percent i)
  (* (/ (width i) (center i)) 100))

練習問題 2.13

全ての数値が正であると仮定する。 { \displaystyle R1 }, { \displaystyle R2 }区間を、{ \displaystyle a }, { \displaystyle c } は中央値, { \displaystyle b }, { \displaystyle d } は誤差とする。区間(下限, 上限) と表現する。また、本文中では「パーセント許容誤差」という単語と「許容誤差」という単語が混在しているが、ここでは許容誤差と統一する。

以下証明。

{ \displaystyle R1 }{ \displaystyle a \pm b (0 \leqq b \leqq 1) } であるとすると、

{ R1 = \left[a(1-b),\ a(1+b)\right] }

{ \displaystyle R2 }{ \displaystyle c \pm d (0 \leqq d \leqq 1) } であるとすると、

{ R2 = \left[c(1-d), c(1+d)\right] }

になる。

この時、練習問題 2.11 で定義した mul-interval区間両方が正である場合の計算から、

{ R1 \cdot R2 = \left[ac(1-b)(1-d), ac(1+b)(1+d)\right] }

になる。

ここで

{ ac(1-b)(1-d) = ac(1 - (b+d) + bd)}

であるが、問題文の仮定から許容誤差 { \displaystyle b }, { \displaystyle d } は非常に小さいため { \displaystyle b \cdot d \fallingdotseq 0} と考えられるので、

{ ac(1 - (b+d) + bd) = ac(1 - (b+d))}

となる。

同様に

{ ac(1+b)(1+d) = ac(1 + (b+d))}

従って区間 { \displaystyle R1 }, { \displaystyle R2 } の積は

{ R1 \cdot R2 = \left[ac(1-(b+d)), ac(1+(b+d))\right]}

ここから中央値は

{ \displaystyle \frac{(ac(1-(b+d)) + ac(1+(b+d)))}{2} = \frac{ac(1+1-(b+d)+(b+d))}{2} = ac}

幅は

{ \displaystyle \frac{(ac(1+(b+d)) - ac(1-(b+d)))}{2} = \frac{ac( 1+(b+d) - (1-(b+d)) )}{2} = ac(b+d)}

最後に、許容誤差は 幅 / 中央値 で求められるので

{ \displaystyle \frac{ac(b+d)}{ac} = (b+d)}

よって二つの区間の積の許容誤差は、因数の許容誤差である { \displaystyle b}{ \displaystyle d} を使って { \displaystyle (b+d)} と表現できることがわかった。

証明終わり。

 実際に計算してみると、以下のようにほぼ (b+d) になる。

;; 1.0 を掛けて実数に変換
gosh> (* (percent (mul-interval (make-center-percent 100 0.00001) (make-center-percent 100 0.00001))) 1.0)
2.0000000004074335e-5

練習問題 2.14

まず適当に par1par2 の結果を見てみる。

gosh> (par1 (make-center-percent 5 1) (make-center-percent 10 5))
(3.0241157556270095 . 3.669550173010381)
gosh> (par2 (make-center-percent 5 1) (make-center-percent 10 5))
(3.2543252595155714 . 3.409967845659164)

gosh> (par1 (make-center-percent 10 10) (make-center-percent 20 15))
(4.5 . 9.730769230769232)
gosh> (par2 (make-center-percent 10 10) (make-center-percent 20 15))
(5.884615384615384 . 7.4411764705882355)

このように全く異なる結果になる。

ある区間 A と B を作成し、式 A/A と A/B の計算においてそれらを用いよ。

(Page 100).

ここで特に A/A に注目してみる。

gosh> (let ([A (make-center-percent 10 10)]
            [B (make-center-percent 20 15)])
        (div-interval A A))
(0.8181818181818182 . 1.222222222222222)

このように、予期していた (1 . 1) から大きくズレている。ここでは、このズレは除算による誤差であると考えた。次に問題文の続きを読むと以下のようにある。

幅が中央値の小さなパーセンテージである区間を用いることで多くの実態を掴むことができるだろう。

(Page 100).

そこで実際に試してみる。

gosh> (let ([A (make-center-percent 10 0.01)])
        (div-interval A A))
(0.9998000199980004 . 1.000200020002)

このようにみると誤差が小さくなったように見える。

 しかし以下の第4回の原稿というものを読んでみるとこの考えは間違っていたことがわかった。

応用数理』チュートリアル 「精度保証付き数値計算」

以下代数学における用語について正確で無い可能性があるが、自分なりに説明を試みてみる。

 上記の原稿の 「2.4 区間演算の性質」において、区間演算では加法と乗法に関する逆元が存在しないとなっている。実際計算してみると、[1,2] - [1,2] = [-1,1] であるが、[-1,1] + [1,2] = [0,3] となり元に戻らないため、加法において [-1,1][1,2]単位元ではない。同様に乗算についても [1,2] / [1,2] = [1/2, 2] となり、これもあきらかに乗法における単位元にはならい。そうなると P99 の式は「代数的」には等価であるのだが、区間である { \displaystyle R_1 }{ \displaystyle R_2 } は逆元を持たないので { \displaystyle \frac{R_1}{R_1} = 1 } とはならない。これが原因で par1par2 の答が一致しないと考えられる。  代数学において、任意の { \displaystyle a \in G } に対して { \displaystyle b \cdot a = a \cdot b = e } を満たす { \displaystyle b \in G } が存在することが群の条件のひとつであるが、区間は加算についても乗算についてもこれを持たない。以上のことから、{ \displaystyle \frac{R_1}{R_1}=1 } といった代数的な演算は適用できないと考えられる。つまり、これは区間が (少なくとも今回の区間の定義においては) 代数的操作をするための性質を持たないことによって起こる問題であると考えられる。

練習問題 2.15

彼女は Alyssa のシステムを用いて区間の計算をする式が、もし式が不確かな値を表現する変数がどれも繰り返されない形であれば、より厳しいエラーの限界を算出すると言う。従って彼女は抵抗の並列に対し、par2 の方が par1 より “より良い” プログラムであると述べた。彼女は正しいだろうか? それは何故か?

(Page 100).

文がいまいちピンと来なかったが、参加者が「不確かな値を表現する変数」というのは並列接続の抵抗の式における { \displaystyle R_1 }{ \displaystyle R_2 } だとすれば、最初の式 (par1) の積 { \displaystyle R_1 R_2 } や、それを { \displaystyle R_1 + R_2 } で除算することのように、変数同士の計算を行うことでこの不確かな変数同士の計算によってより大きな不確かさが生じてしまうということを指しているのではと考えた。これと比べて二番目の式 (par2) では 1 という定数と変数の計算が殆どであり、不確かさが生じそうなのは { \displaystyle \frac{1}{R_1} } の結果と { \displaystyle \frac{1}{R_2} } の結果を足すところのみであるので、こちらの方が不確かさ (=誤差?) が小さくて済みそうである。

 したがって Eva は正しいと言えそうだ。

練習問題 2.16

(警告:この問題はと ても難しい)

(Page 100).

とのことであまり深入りはしていないがぱっと思いついたことを書いておく。

 練習問題 2.14 で言ったように、今回の区間は加法と乗法について逆元を持たないため、通常扱っているような数値と同じような式変形はできない。逆を言えば代数的な変換をするためには逆元を持つように区間というものを表現できれば良いということで、区間というものをそのような形で表現する方法を考えれば良いのだろうか。  練習問題 2.14 で挙げた参考文献にはこれに対する手法が書かれているようだ。なんらかの手法で区間というものを別な次元に写像して、そこでの変換を行うというようなアプローチになるのだろうか。すぐには理解できそうになかったので一旦諦めた。掘り下げたくなったら参考文献などに当たってみたい。