モンテカルロ法で琵琶湖の面積の割合を求める

Disclaimer

このPostには、技術的な内容と同様かそれ以上に滋賀県民が滋賀県をひたすら愛するだけの内容が含まれています。純粋な技術記事を読みたい方は他をあたってください。技術的にすごいことは何も起きていません。

滋賀県民

こんにちは。滋賀県民です。
皆さんは、 1/6という数字を聞いて何を思い浮かべるでしょうか。

実はこれは、滋賀県全体の面積にしめる琵琶湖の面積の割合です。
この数字は、県内の小学校5年生以上なら誰しも必ず知っているはずの値で、結構な割合の県民は知っています。

というのも、滋賀県では、フローティングスクール　という事業が展開されており、県下の小学5年生は、琵琶湖上に浮かぶ学習船(一応県立学校という位置づけです)の中で一泊の宿泊を含む研修活動を行います。この船は、うみのこ　と呼ばれていて、2018年、30年以上親しまれた船が引退して、今はリニューアルされた新船となっています。

なかなかいい体験なので、滋賀県民で良かったことの一つです。

ぜひ一度ご覧ください。 : https://uminoko.jp

ちなみにですが、このうみのこに乗船する小学生が、研修前に覚える「希望の船」という歌を、「地味にまだ歌える」という人が多くて面白いです。無論、私もその一人です。

この船で、県下の小学生たちは、プランクトンの観察やロープワーク体験、カッター(小型の船)での活動などを行います。そのために事前・事後学習を死ぬほどやるのですが、その中でよく出てくる値のうちの一つがこの1/6なのです。

最近では、この 1/6の知名度にあやかって、コクヨが ロクブンノイチ野帳　なるものも出しています。琵琶湖水系沿岸のヨシから生産された紙を使用していて、なかなか良いのでぜひご購入ください。

さて、本題に移りましょう。

ここまでご説明したように、1/6　なのですが、残念ながら京都民や大阪民は口を揃えてこう言います。 「滋賀とか半分くらい琵琶湖ちゃうの?」

そう思うのも当然でしょう。彼らには琵琶湖の周りに広がる豊かで広大な山林が見えていないのでしょう。

しかしながら、流石にそれを甘んじて受け入れるわけにはいきません。

というわけで、大阪、京都民にもできる。というか全世界どこにいてもできる、面積の割合の実験を提案します。

モンテカルロ法

この記事をお読みの方はご存知でしょう。
円周率を求めるのにも非常によく使われる手法です。
円周率を求めるときには、正方形の中にしめる円の面積の割合を求めます。

そう、ある領域における面積の割合は、一定の精度までは確率的に求めることが可能です。
というわけで、滋賀県にしめる琵琶湖の面積をモンテカルロ法で求めます。

地図の入手

まずは、琵琶湖の面積と滋賀県の面積を決定させるベースになる地図を用意します。
とりあえずこれでいいでしょう。

https://www.freemap.jp/itemFreeDlPage.php?b=shiga&s=shiga

地図の処理

次に、白地図に色をつけます。
これは、地図中の 滋賀県ではないエリア 滋賀県の陸のエリア 滋賀県の琵琶湖のエリア をソフトウェアから判別できるようにするための処理です。
適当に塗りつぶしてください。
今回、私はSeashoreを使って 滋賀県の陸のエリア にrgb(255,38,0)、滋賀県の琵琶湖のエリアにrgb(0,0,255)の色をつけました。

ぷよぐらみんぐ

基本的には円周率を求めるのと同じで、領域(正方形とか滋賀県とか)中にランダムな点を設定し、それが指定の領域(円とか琵琶湖とか)に入っているかを判定する処理を繰り返します。
本当に円とやっていることは変わらないので、あまり気にすることないと思いますが、少し気にしないといけないのは、分母はepoch数ではありません。
琵琶湖に入った数を滋賀県に入った数で割ってください。
これを識別するために色をつけています。
画像中のある座標の色のデータは、皆さんご存知のようにOpenCVを用いてimreadして出てくる配列の[y,x]番地にアクセスすることで得られます。
そのほかに気にすることは特にないと思います。
一応これでは試行回数は50000回に設定しています。

import cv2
import numpy as np

print("Biwako Area inspector")
print("©︎ Yusei ito 2020")
print("=========================")


print("Importing image..")
img = cv2.imread('./map.png') #Processed map image file
height, width, channels = img.shape[:3]
print("Image Imported.")
print("\tWidth:"+str(width))
print("\tHeight:"+str(height))

count_max=50000
count_shiga=0
count_biwako=0

for i in range(0,count_max):
    px=img[np.random.randint(height),np.random.randint(width)]#Be attention. Array is formed as [y,x] and return is formed [blue,green,red]
    if (px[0]==0 and px[1]==38 and px[2]==255):
        count_shiga+=1
    if (px[0]==255 and px[1]==0 and px[2]==0):
        count_shiga+=1
        count_biwako+=1
    if (count_max%10==0):
        devisionBy=1 if count_shiga==0 else count_shiga
        print(str(i)+"\t"+str(count_biwako/devisionBy))


print("All Process succeed.\r\n\r\n")    
print("Attempt:"+str(count_max))
print("Shiga:"+str(count_shiga))
print("Biwako:"+str(count_biwako))

print("Overall Score:"+str(count_biwako/count_shiga))

実行結果は、冗長なので一番最後だけ貼ります。
ロクブンノイチなので、期待値は1/6=0.1666666....です。

Attempt:50000
Shiga:21143
Biwako:3493
Overall Score:0.16520834318687036

良いですね。なかなか良いです。

というわけで、滋賀県はいいところです。一度はお越しくださいね。

この投稿に関するもの一式はGitHubにまとめています。合わせてご覧ください。