第5回：■ 条件式・■ 条件分岐

■ 数の大小比較

値が等しい・異なる

演算子 == は、a == b のように用いて、値が等しいか否か判定する。演算子 != は、値が異なるか否か判定する。成立すれば(真ならば) trueが、成立しなければ(偽ならば) false が結果となる。

julia> 1 == 1
true

julia> 1 == 2
false

julia> 1 != 1
false

julia> 1 != 2
true

より大きい

演算子 > は、a > b のように用いて、aの値がbの値よりも大きいか否か。演算子 >= は、a >= b のように用いて、aの値がbの値以上であるか否か判定する。

julia> 2 > 1
true

julia> 2 >= 1
true

julia> 2 >= 2
true

より小さい

演算子 < は、a < b のように用いて、a の値が b の値よりも小さいか否か。演算子 >= は、a <= b のように用いて、a の値が b の値以下であるか否か判定する。

julia> 2 < 1
false

julia> 2 <= 1
false

julia> 2 <= 2
true

■ 論理型

比較演算子の結果は true または false のどちらかである。この二つの値からなるデータの種類を論理型 (logical type)という。

■ if文

Conditional Evaluation

if文は、直後に論理式をとる。論理式の値が true なら、if文の次の文から、end, else, elseif が出現するまでの文を実行する。

条件が成り立つときだけに実行される部分をifブロック(block)という。

ブロックは、字下げ (indent) で表記される。が、字下げは見やすさのためだけである。

julia> x=1
1

julia> y=2
2

julia> if x < y
         println("x は y より小さい")
       end
x は y より小さい

ifブロックの後ろに、else文が続く場合がある。論理式の値が false なら、else文の次の文から end が出現するまでの文 (elseブロック)を実行する。

julia> x=1
1

julia> y=2
2

julia> if x < y
       　println("x は y より小さい")
       else
       　println("x は y より小さくない")
       end
x は y より小さい

else文の前に、elseif文が続く場合もある。最初の if文の論理式が false なら、 elseif文の論理式を計算し、それが true なら、elseif 文の次の文から、elseifまたはend が出現するまでの文 (elseifブロック)を実行する。

julia> x=1
1

julia> y=2
2

julia> if x < y
         println("x は y より小さい")
       elseif x > y
       　println("x は y より大きい")
       else
         println("x は y と等しい")
       end
x は y より小さい

■ if式

if式は、論理式が成立したブロックの最後の値を、式の値とする。

julia> x=40
40

julia> m=if x >= 100
         "x は 100 以上である"
       elseif x >= 50
         "x は 50 以上である"
       elseif x >= 20
         "x は 20 以上である"
       else
         "x は 20 よりも小さい"
       end
"x は 20 以上である"


julia> @show m
m = "x は 20 以上である"
"x は 20 以上である"

if式を用いて ▶︎ 絶対値関数を、以下のように書くこともできる。

julia> myabs1(x) = if x >= 0
         x
         else
         -x
         end
myabs1 (generic function with 1 method)


julia> @show myabs1(-1)
myabs1(-1) = 1
1

■ 3項演算子

Operators

a ? b : c

条件 a が真(true)なら b を、偽なら c を値とする「マクロ」である。

julia> m= 2 > 1 ? "yes" : "no"
"yes"

julia> @show m
m = "yes"
"yes"

if式を用いて ▶︎ 絶対値関数を、以下のように書くこともできる。

julia> myabs2(x) = x >= 0 ? x : -x
myabs2 (generic function with 1 method)


julia> @show myabs2(-1)
myabs2(-1) = 1
1

■ 論理演算

論理否定

Base.:! — Function

論理否定 ! a は、aの論理値を反転する。

julia> ! true
false

julia> ! false
true

論理積

&& - Keyword

論理積 $a\cdot b$は、$a$と$b$との両方が true のときだけ true となる。 $a$と$b$との、少なくとも一つが false なら、false　となる。

julia> true && true
true

julia> true && false
false

julia> false && false
false

julia> false && false
false

論理和

|| - Keyword

論理和論理積 $a + b$は、$a$と$b$との少なくとも一つがtrue のとき true となる。 $a$と$b$とのどちらも false のときは、false　となる。

julia> true || true
true

julia> true || false
true

julia> false || false
false

julia> false || false
false

■ 数の大小比較の連続した記述

数の比較演算子は、連続して記述できる。

julia> 1 < 2 < 3
true

julia> 1 < 2 && 2 < 3
true

▶︎ 総当たりによる不定方程式の解法

方程式の数よりも、未知数の数が少ない方程式を不定方程式という。未知数の性質が決まっていれば、未知数の候補を、方程式に代入して、成り立つ未知数の組合せを求めることができる (総当たり攻撃, brute-force attack)。

問題　「負ではない三つの整数 $x, y, z$ が、次の二つの等式を同時に満たすという。

\[x + y + z = 24\]

\[x+2y+4z = 51\]

三つの数 $x, y, z$ の組合せを全て求めよ。」

ここで、 $x, y, z$ は、0から24までの整数である。 $x, y, z$ の全ての組み合わせに対して、二つの等式が成り立つ条件を、if文に渡す。

julia> for z=0:24
         for y=0:24
           for x=0:24
             if x+y+z==24 && x+2*y+4*z==51
               @show x,y,z
             end
           end
         end
       end
(x, y, z) = (1, 21, 2)
(x, y, z) = (3, 18, 3)
(x, y, z) = (5, 15, 4)
(x, y, z) = (7, 12, 5)
(x, y, z) = (9, 9, 6)
(x, y, z) = (11, 6, 7)
(x, y, z) = (13, 3, 8)
(x, y, z) = (15, 0, 9)

▶︎ 格子点による平面領域の塗り分け

平面座標の第一象限 $0 \le x \le 1, 0 \le y \le 1$ の範囲に、格子点を配置しよう。

そのうち、$x^2+y^2 \lt 1$ の範囲にある点を赤で、それ以外を青で色分けして示そう。

二重の繰り返し (二重ループ)を用いた。

刻み幅は各軸 0.025 に選んだので、格子点の総数は ${40}\times{40}=1600$ 個である。

赤色の範囲は、4分円となる。

using PyPlot
for y=0:0.025:1
  for x=0:0.025:1
    c = x*x + y*y < 1 ? "r" : "b"
    plot(x,y,".", color=c)
  end
end
plt[:axes]()[:set_aspect]("equal")

今度は、二つの不等式 $y < x$ と $y < 1-x$の両方に囲まれる領域を、赤で示そう。

二重ループを一つの for文にまとめた。for文の右側に書かれたループ変数が内側のループに相当する。

using PyPlot
for y=0:0.025:1, x=0:0.025:1
  if y < x && y < 1-x
    c="r"
  else
    c="b"
  end
  plot(x,y,".", color=c)
end
plt[:axes]()[:set_aspect]("equal")

◀︎ 練習

格子点の刻み幅を大きく、または、小さくして、色分けの様子を観察してみよ。

関数plotのfmtパラメータを . から o (circle)に変えてみよ。さらに、markersizeパラメータを加えて、マーカーの大きさを調整できる。

plot(x,y,"o", markersize=3)

Note

刻み幅を小さくすると、格子点の数は急激に増える (例えば、$x,y$軸方向の刻みの数を2倍にすると、格子点の数は 4倍になる)ので、刻み幅を急激に小さくすべきではない。計算量が増えた場合、数分以内に計算が終わらない可能性もあるが、このような状況を経験するのも重要である。

◀︎ 練習

以下の領域を表す条件を、式で表せ。その領域を、格子点を用いて塗り分けてみよ。

点 (0,0)と点(1,0)を一辺とし、第一象限に描かれた正三角形
中心 $\left(\dfrac{1}{2},\dfrac{1}{2}\right)$ , 半径 $\dfrac{1}{2}$ の円

■ 擬似乱数

計算機で発生する乱数を、擬似乱数 (pseudorandom numbers)という。

Note

本来の意味の「乱数」は、計算する方法がないはずである。にもかかわらず、ある手順で乱数を発生しているので「擬似乱数」という。

関数 rand()は、0から1未満の擬似乱数を発生する。

julia> rand()
0.7464052231559362

julia> rand()
0.818784058672192

julia> rand()
0.9689823232186385

Base.Random.srand - Function

関数 srand(m)は、擬似乱数の種をリセットする。種mには、0から$2^{32}-1$までの整数を指定する。同じ種を指定すると、同じ系列で擬似乱数を発生する。

julia> srand(1234)
MersenneTwister(UInt32[0x000004d2], Base.dSFMT.DSFMT_state(Int32[-1393240018, 1073611148, 45497681, 1072875908, 436273599, 1073674613, -2043716458, 1073445557, -254908435, 1072827086  …  -599655111, 1073144102, 367655457, 1072985259, -1278750689, 1018350124, -597141475, 249849711, 382, 0]), [1.38663, 1.2279, 1.59552, 1.56672, 1.35442, 1.87451, 1.2136, 1.50203, 1.49049, 1.4059  …  1.70015, 1.56379, 1.60416, 1.46075, 1.86963, 1.6698, 1.70239, 1.44344, 1.927, 1.92703], 382)

julia> rand()
0.5908446386657102

julia> rand()
0.7667970365022592

julia> rand()
0.5662374165061859

▶︎ 乱数による平面領域の塗り分け

関数 rand() を2回用いて、座標点 (x,y) を発生しよう。点の数を増やすと、これらの点は、平面の第一象限 $0 \le x \lt 1, 0 \le y \lt 0$ の範囲を埋め尽くすはずである。

そのうち、$x^2+y^2 \lt 1$ の範囲にある点を赤で、それ以外を青で色分けして示そう。

点の個数は $2^{10}=1024$ である。格子点で図示した場合 ( ▶︎ 格子点による平面領域の塗り分け ) に比べて点の数は少ないが、特徴を捉えている。

using PyPlot
n=2^10
for i=1:n
  x=rand()
  y=rand()
  c = x*x + y*y < 1 ? "r" : "b"
  plot(x,y,".", color=c)
end
plt[:axes]()[:set_aspect]("equal")

今度は、二つの不等式 $y < x$ と $y < 1-x$ の両方に囲まれる領域だけ、赤で示そう。

using PyPlot
n=2^10
for i=1:n
  x=rand()
  y=rand()
  if y < x && y < 1-x
    c="r"
  else
    c="b"
  end
  plot(x,y,".", color=c)
end
plt[:axes]()[:set_aspect]("equal")

▶︎ モンテカルロ法による平面図形の面積の推定

以上の例で、ランダムに落とした点の総数のうち、図形の中に入った数を数えよう。両者の割合から、図形の面積を推定できる。これを、モンテカルロ (Monte Carlo)法による面積の算出法という。

Note

モンテカルロ (Monte Carlo)は、モナコ公国(Monaco)の地名の一つであり、公営カジノで有名である。モンテカルロ法 (Monte Carlo Method)は、擬似乱数を用いて、数値計算やシミュレーションを行う手法の総称である。(面積を推定する手法以外の「モンテカルロ法」もある)

まず、4分円の面積を推定しよう。

n=2^10
s=0
for i=1:n
  x=rand()
  y=rand()
  if x*x + y*y < 1
    s += 1
  end
end

s / n = 0.7685546875
pi / 4 = 0.7853981633974483
(s / n - pi / 4) / (pi / 4) = -0.021445779583424734
-0.021445779583424734

4分円の面積は $\dfrac{\pi}{4}$である。点数 1024個で、相対誤差 2% 程度の、面積推定値が得られた。

今度は、点の総数に対して、面積の推定値を描いてみる。

using PyPlot
for m in 1:16
  n=2^m
  s=0
  for i=1:n
    x=rand()
    y=rand()
    if x*x + y*y < 1
      s += 1
    end
  end
  plot(n, s/n, ".")
  @show n, s/n
end
ylim(0.9*pi/4, 1.1*pi/4)
xlabel("n")
xscale("log")
axhline(pi/4, color="k", lw=0.5)

(n, s / n) = (2, 0.5)
(n, s / n) = (4, 1.0)
(n, s / n) = (8, 1.0)
(n, s / n) = (16, 0.625)
(n, s / n) = (32, 0.84375)
(n, s / n) = (64, 0.75)
(n, s / n) = (128, 0.796875)
(n, s / n) = (256, 0.78125)
(n, s / n) = (512, 0.76953125)
(n, s / n) = (1024, 0.7900390625)
(n, s / n) = (2048, 0.77734375)
(n, s / n) = (4096, 0.7783203125)
(n, s / n) = (8192, 0.7891845703125)
(n, s / n) = (16384, 0.78314208984375)
(n, s / n) = (32768, 0.79083251953125)
(n, s / n) = (65536, 0.7855987548828125)

正しい面積との相対誤差を、点の総数に対して描いてみる。

using PyPlot
for m in 1:16
  n=2^m
  s=0
  for i=1:n
    x=rand()
    y=rand()
    if x*x + y*y < 1
      s += 1
    end
  end
  plot(n, abs(s/n-pi/4)/(pi/4), ".")
end
xlabel("n")
ylabel("relative errors")
xscale("log")
yscale("log")

◀︎ 練習

モンテカルロ法を用いて、以下の領域の面積を推定せよ。(これらの領域の、格子点を用いた塗り分けは、以前の例題・練習として行った)。

余裕があれば、点数に対する相対誤差の変化も描いてみよ。

二つの不等式 $y < x$ と $y < 1-x$ の両方に囲まれる領域
点 (0,0)と点(1,0)を一辺とし、第一象限に描かれた正三角形
中心 $\left(\dfrac{1}{2},\dfrac{1}{2}\right)$ ,　半径 $\dfrac{1}{2}$ の円

■ 関数の定義 (代入文形式)

関数は、いくつかの値を受け取って、何らかの操作をして返す、まとまった処理である。

1行によるユーザ定義関数の例を、以下に示す。

julia> f(x)=x*2-1
f (generic function with 1 method)

(の前が、ユーザ定義関数の名前である。関数名の規則は、■ 変数名の規則と同じである。

括弧 () の中に、変数名の並びを記述する。この変数名を、仮引数(「かり・ひきすう」 parameter)という。

=より右が、関数の定義である。具体的な値（実引数「じつ・ひきすう」 argument)を入れて、関数を評価する (evaluate)。関数の定義の中に出現した仮引数は、実引数の値に置き換えられる。

julia> f(0)
-1

julia> f(1)
1

実引数にベクトルなどを与えて、各要素に対して評価する場合には、関数名の直後にピリオド .を置く。

julia> f.([1,2,3])
3-element Array{Int64,1}:
 1
 3
 5

julia> f.(0:5)
6-element Array{Int64,1}:
 -1
  1
  3
  5
  7
  9

仮引数は、2個以上でもよい。

julia> g(x,y)=x*y
g (generic function with 1 method)

julia> g(1,1)
1

julia> g("a "," b")
"a  b"

数同士の * 演算は、数の乗算である。文字列同士の *演算は、文字列の連結である。

引数の型にあわせて、正しい演算が得られていることに着目せよ。

▶︎ 「はさみうち」法による、方程式の求解

条件判断を、繰り返し行うことで、求める答えに近づいて行く例を、もう一つ紹介する。

▶︎ 不連続な有理式を描くの例として挙げた分母の式 $f(x)=x^3+3x^2-4x-12$ について、方程式 $f(x)=0$ の解の近似値を求めてみる。

まず、$y=f(x)$のグラフを描こう。

f(x)=x^3+3x^2-4*x-12
xs=-3.5:0.05:3
#
using PyPlot
plot(xs, f.(xs))
axhline(0, color="k", lw=0.5)

さて、連続な関数 $f(x)$ が、区間 $a \lt x \lt b$ で単調（単調増加または単調減少）であるとする。ここで、$f(a)$ と $f(b)$ との符号が異なるとき $f(a)\cdot f(b)\lt 0$ 、方程式 $f(x)=0$ の解は $a \lt x \lt b$ の区間にある。

ここで、$a$ と $b$ との中点 $c=\dfrac{a+b}{2}$ をとり、$f(a)$ と $f(c)$ が同じ符号であれば左端$a$を$c$に更新する。逆に、$f(b)$ と $f(c)$ が同じ符号であれば右端$b$を$c$に更新する。

この手順を繰り返すことで、方程式 $f(x)=0$ の解が存在する区間 $a \lt x \lt b$ を狭めていくことができる。この手法を「はさみうち」法という。英語では、squeeze theorem, pinching theorem, sandwich theorem などと呼ばれる。

上の関数 $f(x)$ で $a=-3.2 < x < b=-2.6$ の区間を選ぼう。

この区間で $f(x)$ が単調であることを、まず確かめよう。

# f(x)=x^3+3*x^2-4*x-12
xs=-3.2:0.01:-2.6
plot(xs, f.(xs))
axhline(0, color="k", lw=0.5)

この区間で、「はさみうち」を数回繰り返してみる。

# f(x)=x^3+3x^2-4*x-12
a=-3.2; b=-2.6
@show a,b
@show f(a),f(b)
for i=1:10
  c=(a+b)/2
  @show i, a, b, c, f(c)
  if f(a)*f(c) > 0
    a = c
  else
    b = c
  end
end

(a, b) = (-3.2, -2.6)
(f(a), f(b)) = (-1.248000000000001, 1.1039999999999974)
(i, a, b, c, f(c)) = (1, -3.2, -2.6, -2.9000000000000004, 0.4409999999999954)
(i, a, b, c, f(c)) = (2, -3.2, -2.9000000000000004, -3.0500000000000003, -0.26512500000000294)
(i, a, b, c, f(c)) = (3, -3.0500000000000003, -2.9000000000000004, -2.9750000000000005, 0.12126562499999949)
(i, a, b, c, f(c)) = (4, -3.0500000000000003, -2.9750000000000005, -3.0125, -0.06343945312499955)
(i, a, b, c, f(c)) = (5, -3.0125, -2.9750000000000005, -2.9937500000000004, 0.031015869140624375)
(i, a, b, c, f(c)) = (6, -3.0125, -2.9937500000000004, -3.0031250000000003, -0.015683624267579077)
(i, a, b, c, f(c)) = (7, -3.0031250000000003, -2.9937500000000004, -2.9984375000000005, 0.007797855377194907)
(i, a, b, c, f(c)) = (8, -3.0031250000000003, -2.9984375000000005, -3.00078125, -0.003909912586216535)
(i, a, b, c, f(c)) = (9, -3.00078125, -2.9984375000000005, -2.9996093750000004, 0.0019522095322592747)
(i, a, b, c, f(c)) = (10, -3.00078125, -2.9996093750000004, -3.0001953125000003, -0.000976791389289744)

この範囲の解は $x=-3$である。

@show f(-3)

f(-3) = 0
0

区間が狭まる様子を、グラフに描いてみよう。

# f(x)=x^3+3x^2-4*x-12
using PyPlot
a=-3.2; b=-2.6
for i=1:15
  plot(b,i, "ro")
  plot(a,i, "b.")

  c=(a+b)/2
  if f(a)*f(c) > 0
    a = c
  else
    b = c
  end
end
xlabel("x")
ylabel("i")

今度は、区間の幅 $b-a$をプロットしよう。区間の幅が、単調に減少する様子が観察された。

# f(x)=x^3+3x^2-4*x-12
using PyPlot
a=-3.2; b=-2.65
for i=1:30
  plot(i, b-a, "r.")

  c=(a+b)/2
  if f(a)*f(c) > 0
    a = c
  else
    b = c
  end
end
yscale("log")
ylabel("b-a")
xlabel("i")

◀︎ 練習

同じ関数で、異なる区間で、はさみうち法を試してみよ。
別の関数に対して、単調な区間を選び、はさみうち法を試してみよ。

▼◀︎ NLsolveパッケージの紹介

NLsolveパッケージは、連立非線形方程式を解くためのパッケージのひとつである。

NLsolveパッケージを用いて、上の例題を解く例を示す。ただし、まだ説明していない文法(ベクトルの添字、関数の副作用)を用いているので、学習が進んでから、再度試してみよ。

julia> # Pkg.add("NLsolve") # パッケージの導入。冒頭の`#`を外して１回だけ実行すればよい。
       
       using NLsolve


julia> # 解くべき関数を定義する
       function g!(r, x)
           r[1] =x[1]^3+3*x[1]^2-4*x[1]-12
       end
g! (generic function with 1 method)


julia> # nlsolve関数を呼び出す
       nlsolve( g!, [-3.5], autodiff = :forward)
Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [-3.5]
 * Zero: [-3.0]
 * Inf-norm of residuals: 0.000000
 * Iterations: 5
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 6
 * Jacobian Calls (df/dx): 6

★ 今回のまとめ

数の大小比較
if文
if式
3項演算子
論理演算
格子点による平面領域の塗り分け
擬似乱数
乱数による平面領域の塗り分け
モンテカルロ法による平面図形の面積の推定
ユーザ定義関数 (代入文の形式)
「はさみうち」法による、方程式の求解