MCMCへの道（２）

 MCMCへの道（１）からの続き

さて，MCMC法#2棄却サンプリングである。手順は次の通り。

(1) 目標分布 p(x) を事後分布ともいう。これは既知の関数。

(2) 乱数発生計算が容易な分布 q(x) で， k q(x) ≧ p(x) を満たすものを選ぶ。 kは正定数。例えば，一様分布とか正規分布。

(3) q(x) にしたがう乱数 x' を発生する。

(4) [0, k q(x)] の一様分布にしたがう乱数 y' を発生する。

(5) y' ≦ p(x') ならx' を採用する。そうでなければ棄却する。

この(3) から(5) を必要なだけ繰り返す。

using Plots

function beta(x,a,b)
  return x^(a-1)*(1-x)^(b-1)
end

x= [k for k in 0.0:0.001:1.0]
y= beta.(x,10.2,5.8)
z= [0.00013 for k in 0.0:0.001:1.0]
scatter(x,y,xlim=(0,1.0),ylim=(0,0.00015),
  markershape = :circle,
  markersize = 1,
  markeralpha = 0.75,
  markercolor = :blue,
  markerstrokewidth = 0.1,
  markerstrokealpha = 0.1,
  markerstrokecolor = :white)

g1=scatter!(x,z,xlim=(0,1.0),ylim=(0,0.00015),
  markershape = :circle,
  markersize = 1,
  markeralpha = 0.75,
  markercolor = :orange,
  markerstrokewidth = 0.1,
  markerstrokealpha = 0.1,
  markerstrokecolor = :white)

n=3000
p=rand(n)
r=0.00013*rand(n)
q=ones(n)
for i in 1:n
  f=beta(p[i],10.2,5.8)
  q[i]=ifelse(f <= r[i],0,r[i])
end

g2=scatter(p,q,xlim=(0,1.0),ylim=(0,0.00015),
  markershape = :circle,
  markersize = 1,
  markeralpha = 0.75,
  markercolor = :green,
  markerstrokewidth = 0.1,
  markerstrokealpha = 0.1,
  markerstrokecolor = :white)

plot(g1,g2,legend=:left,aspect_ratio=6000)

図：ベータ分布 p(x) と一様分布 q(x) と定数 k=0.00013 による棄却法の例

女紅場

 法事で京都の街を歩いていたら，丸太町の鴨川縁に女紅場（にょこうば）の碑があった。京都府立第一高女がなんとかとあったので，調べてみた。

女紅場とよばれた女子の教育機関にはいくつかの類型があるようだが，この京都の女紅場は1872年に設置された日本最初の公立の女学校であり，後の京都府立第一高等女学校の前身だった。現在は京都府立鴨沂高等学校である。

大阪でこれに対応するのが，1874年に堺県に設置された女紅場であり，後の大阪府立泉陽高等学校だ。ここは大阪で２番めに古い高等女学校となった。なお，大阪で最も古い女学校というのは，1882年に設置された府立大阪師範学校の附属裁縫場を起源として大阪府女学校となった，現在の大阪府立大手前高等学校である。

写真：丸太町鴨川西岸にある女紅場の碑（Wikipediaより引用）

［１］女学校（Wikipedia） 

定偏角分散プリズム

ペラン・ブロッカなのかペリン・ブロッカなのかペロン・ブロッカなのかよくわからないが，定偏角分散プリズムというものがあることがわかった。いやあまりわかっていない。幾何光学をまともに勉強していないので。

とりあえず，イメージを掴むために，Mathematicaでプリズムの光路シミュレーションプログラムを書いてみた。わかったような，わからないような・・・それでも高等学校で学ぶ幾何光学の法則さえ知っていればコードは書ける。

プリズムは四角形であり原点(0,0)の内角は90度（直角），y軸上の点（0,2√3-2）は75度，(√3,1)にある点は135度，x軸上の点（√3+1/√3,0）は60度となっていて，その屈折率をnとする。

aは入射光の入射面からはかった入射角度，tbは入射光が屈折した光線の進行方向のプリズム底面に対する傾き，(x1, y1)は全反射する点の座標，tcは反射光の底面に対する傾きの絶対値，tdはプリズムからの出射光の底面からはかった出射角度の勾配，x2は出射点の座標である。

n = 1.6; a = 36.8699;
t1 = Tan[15 Degree]; t2 = Tan[30 Degree]; x3 = Sqrt[3] + 1/Sqrt[3];
tb[n_, a_] := Cos[a Degree]/Sqrt[n^2 - Cos[a Degree]^2]
x1[n_, a_] := Sqrt[3]*t1/(t1 + tb[n, a])
y1[n_, a_] := 1 + Sqrt[3]*t1*tb[n, a]/(t1 + tb[n, a])
tc[n_, a_] := (t2 + tb[n, a])/(1 - t2*tb[n, a])
td[n_, a_] := Sqrt[tc[n, a]^2 - n^2 + 1]/n
x2[n_, a_] := x1[n, a] + y1[n, a]/tc[n, a]

ro = Point[{(1 + 2*t1)/(1 + t1), (1 + 2*t1)/(1 + t1)}];
pol = Polygon[{{0, 0}, {0, 1 + t1}, {1, 1}, {x3, 0}}];

lin[n_, a_] := Line[{{- Sin[a Degree], 1 - Cos[a Degree]}, {0, 1.0}, {x1[n, a], y1[n, a]}, {x2[n, a], 0}, {x3, -td[n, a]*(x3 - x2[n, a])}}]
das[n_, a_] := Line[{{x2[n, a], 0}, {-0.5, -td[n, a] (-0.5 - x2[n, a])}}]

f[n_, a_] := Graphics[{LightGray, pol, Thick, Red, lin[n, a], Dashed, Blue, das[n, a],PointSize[0.03],ro}]

{f[n, a], td[n, a]/Tan[a Degree]}
Manipulate[{f[n, a], n, a, td[n, a]/Tan[a Degree]}, {a, 25, 50}]

図　Pellin-Broca PrismのMathematicaシミュレーション

なお，td[n,a]/Tan[a Degree] = 1ならば入射光と出射光が直交することになる。

［１］定偏角分散プリズムの最適回転中心の位置（小穴純）

人間国宝

例年このころに人間国宝（文化財保護法の重要無形文化財の各個認定保持者）の新規認定が発表される。今年は，人形浄瑠璃文楽人形の桐竹勘十郎さんが入っていた。

Wikipediaで調べてみると，文楽ではこれまでに次の方々が人間国宝に認定されていた。なお，初回は1955年である。文楽協会のページでは，太夫：三味線：人形＝21：21：43で人形遣いの人数が多いのだけれど，人間国宝になるのは難しい。（以下は生年-人間国宝認定年-没年 を表す）

人形浄瑠璃文楽太夫
1 豊竹山城少掾（1878-1955-1967）
2 八世竹本綱太夫（1904-1955-1969）
3 六世竹本住太夫（1886-1955-1959）
4 十世豊竹若太夫（1888-1962-1967）
5 四世竹本津太夫（1916-1973-1987）
6 四世竹本越路太夫（1914-1971-2002）
7 七世竹本住太夫（1924-1989-2018）
8 九世竹本綱太夫→九世竹本源太夫（1932-2011-2015）
9 八世豊竹嶋太夫（1932-2015-2020）
10 豊竹咲太夫（1944-2019-）

人形浄瑠璃文楽三味線
1 四世鶴澤清六（1889-1955-1960）
2 六世鶴澤寛治（1887-1962-1974）
3 野澤松之輔（1902-1972-1975）
4 二世野澤喜左衛門（1891-1962-1976）
5 十世竹澤彌七（1910-1972-1976）
6 四世野澤錦糸（1917-1988-1988）
7 五世鶴澤燕三（1914-1985-2001）
8 八世竹澤團六→七世鶴澤寛治（1928-1997-2018）
9 鶴澤清治（1945-2007-）

人形浄瑠璃文楽人形
1 二世桐竹紋十郎（1900-1965-1970）
2 二世桐竹勘十郎（1920-1982-1986）
3 初代吉田玉男（1919-1977-2006）
4 吉田文雀（1928-1994-2016）
5 吉田簑助（1933-1994-）
6 吉田和生（1947-2017-）
7 三世桐竹勘十郎（1953-2021-）

なお，人形浄瑠璃文楽座は，重要無形文化財の芸能分野の総合認定を受けており，ユネスコの無形文化遺産にもなっている。大阪府・市の維新行政からはひどい扱いを受けているが・・・

カレル大学

石橋や池田に住んでいたころ，あるいは天王寺に通っていたころは，仕事帰りに本屋に寄るのが日課だった。天理に引っ越してからも，大和八木にはかろうじて本屋らしきものがあったので，毎日新刊をチェックすることができた。

でも本屋の消滅が進み始めたころに定年となったたため，本屋にいくのは２，３ヶ月に１度のペースになってしまった。アマゾンやグーグルでは目的の本を探せるかもしれないが，目的としない本に出会うことができない。ネット上でこれを可能にするにはVRが次の段階まで進化して，アーサー・C・クラーク（1917-2008）のダイアスパーの世界が実現するのを待たなければならない。

そんなとき，YouTubeの読書系チャンネルが未知の本を探すための補完になるかもしれない。ヨビノリたくみと斎藤明里のほんタメをみていたら，【全10冊】小説1000冊読んだ私が最近読んだ本というタイトルで，アンナ・ツィマ（1991-）のシブヤで目覚めてを先月末に紹介していた。そうか，深緑野分によればコニー・ウィリス（1945-）なのか。ますます読んでみたくなる。

そのアンナ・ツィマと訳者の阿部賢一・須藤輝彦を招き，実践女子大学文学部准教授のブルナ・ルカーシュが司会をする特別トーク・イベントがあった。アンナ・ツィマとブルナ・ルカーシュが卒業しているのが，プラハにあるカレル大学だ。1348年創設のドイツ語圏最古の大学である。

カレル大学は知らなかったけれど，ヤン・フス，ニコラ・テスラ，カレル・チャペック，フランツ・カフカ，ミラン・クンデラや，多くのノーベル賞学者を輩出しているヨーロッパの主要大学のひとつだった。アルバート・アインシュタインやエルンスト・マッハもこの大学の教員だった。

写真：歴史地区にある旧校舎＝チェコ国立図書館（Wikipediaから引用）

MCMCへの道（１）

 そういえば，機械学習をまともに勉強していなかった。いや，そもそも統計学も確率論もあれもこれもである。先日の水素原子波動関数の可視化についても進んでいないのだけれど，どうやら，メトロポリス・ヘイスティング法を使うのが望ましいらしい（by tsujimotter）。そういえば，モンテカルロ法もまともに勉強してこなかった。

そこで，反省してMCMC（マルコフ連鎖モンテカルロ法）の全体像をマスターすべく，いろいろ探したが，おこちゃま向けの適当なテキストがない。こうなるとYouTubeだのみだ。機械学習基礎理論独習というサイトがあったので，このPythonコードをJuliaで再現しながら学ぶことにした。

その１回目はMCMC法#1モンテカルロ法なので，これを再現してみた。

using BenchmarkTools
using Random
using Plots

function findpi(n)
  rng = MersenneTwister(0)
  count = 0
  for i in 1:n
    count += ifelse(rand(rng)^2+rand(rng)^2<=1.0,1,0)
  end
  return 4*count/n
end

#@btime findpi(10^6)
#@benchmark findpi(10^8)

function circle(n)
  x=rand(n)
  y=rand(n)
  c=zeros(n)
  for i in 1:n
    c[i]=ifelse(x[i]^2+y[i]^2<=1.0,0.0,1.0)
  end
  scatter(x,y,
  marker_z = c,
  markershape = :circle,
  markersize = 1,
  markeralpha = 0.2,
  markercolor = :grey,
  markerstrokewidth = 0.1,
  markerstrokealpha = 0.1,
  markerstrokecolor = :white,
  markerstrokestyle = :dot,
  aspect_ratio = 1.0,
  xlim=(0.0,1.0),
  ylim=(0.0,1.0),
  legend=:topright
  )
end

circle(10000)

jupyterlabを使っているが，多数のデータの描画で表示が重くなってかたまる現象がみられる。黒木さんは400msで10^8点のデータが計算できるとしていたが，こちらではメルセンヌツイスターアルゴリズムを用いて300msくらいかかってしまう。

図：単純なモンテカルロ法による円周率の計算

知の巨人

先日，立花隆（1940-2021）がなくなったときに， NHKをはじめてとしてマスコミが一斉に知の巨人というフレーズで紹介していたのが気持ち悪くて仕方なかった。たぶん田中角栄からはじまる政治・社会問題で名を上げた後に，科学の領域の評論に手を広げたことで，文系のジャーナリストたちが一目置いた的なことが発端のような気がする。あほらしい。

立花隆は，たしかにジャーナリストの嗅覚を働かせて様々な分野に興味を持ったのだけど，あくまでも表層をなぞっていただけとみえる。著書やNHK特集のタイトルとあらすじをみているとわかる。そもそも田中角栄の案件でさえ，アメリカの日本に対する謀略の中の一つの役割を担ったものであり，それが日本にとってよかったのかどうだか。

それでは，知の巨人にふさわしいのは誰かと考えてみたが，原子爆弾や電子計算機の開発や数学・量子力学の基礎論に重要な寄与をしたジョン・フォン・ノイマン（1903-1957）くらいであまり思い浮かばなかった。たまたま，田中康夫の横浜市長選挙出馬にかかわって茂木健一郎が話しているのをきいて，哲学者・論理学者・数学者・社会批評家・政治活動家のバートランド・ラッセル（1872-1970）を忘れていたことに気づいた。

高等学校では理数科だったため，社会の時間は少なくて，1年地理，2年世界史，3年倫理という選択をした。3年のときの倫理の先生は眼鏡の細川先生だった。1971年なのでちょうどバートランド・ラッセルが亡くなったところだったが，授業に分厚い箱入りのみすず書房の西洋哲学史：市井三郎（1922-1989）訳の本を持ってきていて，なんだかすごいと感動したことがある。

細川先生は，学校を抜け出してとなりの食堂から帰ってきた生徒たちを発見しても，そのまま沈黙して通過させてしまうようなことで，寡黙で特徴のない方だったが，倫理の授業というか教科書は非常におもしろくて，現代国語の次に集中して受講した科目だった。

写真：バートランド・ラッセル（ノーベル文学賞のページから引用）

［１］Bertrand Russel（Stanford Encyclopedia of Philosophy）

三島喜美代

 晩夏・小暑・蓮始開（はすはじめてひらく）

先週の日曜美術館で，大阪の十三に住んでいる三島喜美代が取り上げられていた。88歳という年齢を感じさせない意欲を持つパワフルなおばちゃんだった。日本語のWikipediaページはないが，英語版はあるので訳してみた。

三島喜美代（1932-）は、日本の現代美術家であり，新聞や漫画、箱などの「壊れやすい印刷物」を陶土で再現した作品で知られている。1960年代初頭に画家として活動を開始した三島は、1971年から陶芸を始めた。 この頃から，新聞や広告ポスターのイメージをシルクスクリーンの技法で粘土に刷り込むようになる。製造物を用いた作品は、クレス・オルデンバーグやアンディ・ウォーホルの作品、戦後の日本の具体や独逸美術協会の作品と類似している。

三島喜美代は1932年に大阪で生まれ，1951年に大阪の公立扇町高校を卒業した。1986/87年、ロックフェラー奨学金ACCの助成を受けてニューヨークに留学。現在は大阪に在住。

夫の三島茂司（1920-1985）も画家であり，二人展が最近開かれていた。岐阜県の土岐市のアトリエはほとんどゴミの集積場のようだ。

東京の倉庫街にあるART FACTORY城南島に展示されている陶器製の新聞の束は天井まで高く積み上げられて迷路になっている。また，1万個の新聞記事を印刷したレンガが敷き詰められる作品を見て，１つ3000円でも3千万円かと思ったが，市場価格は1つ30万円のオーダーだった。

この情報の化石は，ハードディスクや半導体と違って，長く残り続けるのかもしれない。

写真：三島喜美代 《作品 21-A》2021年

陶にシルクスクリーン印刷，鉄（婦人画報から引用）

コラッツ予想

 コラッツ予想に1億2千万円の懸賞金がかかった。胴元は音圧爆上げくんという音楽系ウェブサービス（料金が無料！）の会社だとのこと。クレイ数学研究所の７つのミレニアム懸賞問題（１件100万ドル）より懸賞金が高い。

それはいいとして，コラッツ予想とは次のようなものだ。ある自然数 n (n≧2) から出発してその次の自然数を，n が偶数のときは n/2 ，n が奇数のときは 3n+1 で定義する。こうして得られた自然数の列が最終的には1に到達するという予想である。

まだ証明はできていないが（だから懸賞金問題になるわけで），2^68 ≒ 3x10^20 まではコンピュータで確かめられている。奥が深そうなので，詳細には立ち入らずに，Juliaでコーディングしてみた。まったく工夫のないプログラムによって，1〜10^8までの奇数に対して1分以内で計算することができた。

なお，プログラミング教育と称して小中学校でモソモソする内容よりは，エクセルを使えるようにするほうが圧倒的に価値があるという説にはほとんど賛同したくなる今日このごろである。

using Plots

function collatz(m,n,p)
# m : number of call
# n : argument of function
# p : print switch
#
  if p==1
    print(n,"->")
  end
  m=m+1
  if mod(n,2)==0
    return m,div(n,2)
  else
    return m,3*n+1
  end
end

function sequence(m,n,p)
# n0: starting argument
# m : iteration number
#
  n0=n
  while n!=1
    m,n=collatz(m,n,p)
  end
  if p==1
    println(1,":",m)
  else
    return n0,m
  end
end

function counter(N,M)
  for i in 1:2:N
    n,m=sequence(0,i,0)
    if m > M
      println(n,":",m)
    end
  end
end

function plotter(N)
  x,y = zeros(N),zeros(N)
  for i in 1:2:N
    x[i],y[i]=sequence(0,i,0)
  end
  scatter(x,y,legend=:topleft)
end

@time counter(100000000,800)
plotter(1000)

図　コラッツ予想 奇数 n=3〜999までのシーケンス数の散布図

一読・十笑・百吸・千字・万歩

午前の昼寝から起きると， 高齢者の健康法というチラシが眼前に届けられました。杏林大学名誉教授石川恭三先生ご指導による日本医師会の健康プラザNo.481（日医ニュース）です。

一読：一日に一度はまとまった文章を読もう。✓（これはクリアできているかな，でもまとまった文章ってどのくらいの量をさすのだろうか，どうやら新聞の社説程度ということで，日経だと2000字弱だ。私の履歴書も同程度なので，文化欄読破というのがよさそうだ）

十笑：一日に十回くらいは笑おう。✓（家族LINEをみていればクリアできそうである）

百吸：一日に百回くらい（一度には十回くらい），深呼吸しよう。×（これはできていません）

千字：一日に千字くらいは文字を書こう。？（認知機能を高めるためにできるだけ漢字を使えとあるので，キーボードはカウントしないのかもしれない。なお千字というのは結構多くてこの日々のブログの平均値ではちょっと足りない）

万歩：一日に一万歩を目指して歩こう。△（目指してはいるのだけれど，毎日の朝の散歩は年平均で4500歩である。歩く日は45分以上で5000歩を超えるのだけれど雨の日もあるからなあ・・・）

図：けんぽれん（健康保険連合会組合）のページから引用

球面調和関数（４）

 球面調和関数（３）からの続き

多くの疑問はよく探してみると答えになるウェブサイトがみつかるものだ。Juliaによる球面調和関数というか原子波動関数の可視化のための関数定義もここにある。面倒なのがラゲール陪関数の定義や規格化なのだけれど，とりあえずなんとかなった。あとは球面調和関数で磁気量子数が負の場合の扱いなど。

using GSL
using Plots
using LinearAlgebra

function Y(l,m,θ,ϕ)
# spherical harmonics
# l: orbital qn, m: magnetic qn
# θ,φ: spherical angular coordinate
#
  am=abs(m)
  if m <= 0 || mod(am,2) == 0
    p=1
  else
    p=-1
  end
#
  if am > l
    return 0
  else
    return p*sf_legendre_sphPlm(l,am,cos(θ))*exp(im*am*ϕ)
  end
end

function R(n,l,r)
# Z: nuclear charge (set 1 for hydrogen)
# n: principle qn, l: orbital qn
# r: Bohr radius unit a0 = ℏc/α・mc^2
#
Z=1
ρ=2*Z*r/n
  if n < l
    return 0
  else
    return sf_laguerre_n(n-l-1,2*l+1,ρ)*exp(-ρ/2)*ρ^l
  end
end

function norm(n,l)
# square of radial wave function norm
  return (2/n)^3 *factorial(n-l-1)/(2*n*factorial(n+l))
end

x=(0.0:0.1:25)
#
r10=(R.(1,0,x).*x).^2 .*norm(1,0)
r20=(R.(2,0,x).*x).^2 .*norm(2,0)
r21=(R.(2,1,x).*x).^2 .*norm(2,1)
r30=(R.(3,0,x).*x).^2 .*norm(3,0)
r31=(R.(3,1,x).*x).^2 .*norm(3,1)
r32=(R.(3,2,x).*x).^2 .*norm(3,2)
#
plot(x,r10,linewidth=1,legend=:topright)
plot!(x,[r20,r21],linewidth=1.5)
g1=plot!(x,[r30,r31,r32],linewidth=2,linestyle=:dash)

t1=(0.0:0.01:pi)
t2=(pi:0.01:2pi)
#
y10=abs2.(Y.(1,0,t1,0))
y11=abs2.(Y.(1,1,t1,0))
y20=abs2.(Y.(2,0,t2,0))
y21=abs2.(Y.(2,1,t2,0))
y22=abs2.(Y.(2,2,t2,0))
#
plot(t1,[y10,y11],proj=:polar,linewidth=2)
g2=plot!(t2,[y20,y21,y22],proj=:polar,linewidth=2)

plot(g1,g2)

図　水素原子の動径部分と角度部分

ワクチン接種統計のマジック

オリンピック開始まであと残りわずかとなった。今日，緊急事態宣言が東京・沖縄で8/22まで発出され，蔓延防止等重点措置が埼玉・千葉・神奈川・大阪で 8/22まで延長された。また，北海道・愛知・京都・兵庫・福岡については7/11までで蔓延防止等重点措置が解除された。

7時のNHKニュースの枠で，菅首相記者会見がだらだらと放送されている。気になったことがある。ワクチン接種のロジスティックスというか管理態勢は破綻していると思うのだが，総数としては必要分が確保されていて（地方自治体に滞留している）というところまではよい。ところで，現時点で1日あたり130-140万回のペースで接種が進んでいると菅は発言した。

いや，いくらなんでもそれはないだろうと，首相官邸のワクチンのページをみると，なにやら内閣官房の新型ワクチンの統計のページの形式が変更されている。あやしい。しかも，新しく設けられている日別のファイルをみると（データ列が途中で入れ替わっているという単純ミスがある），確かに日平均130万回以上接種がすすんでいるかのように見える。

で，よく調べてみると，これには，毎日の遡及入力が含まれない段階の各報告日の過小評価した合計の一覧がならんでいる。一方，それらの最終日において一定の遡及入力が加わったものを並べた結果とは違ってくる。このため，政府発表の数値は，最終的に遡及入力が収束するものと比べて，日々の増加分を過大評価することになっている。グラフでみたほうがわかりやすいだろう。

図　ワクチン接種統計のマジック（2021.7.8の政府統計より）

青が遡及入力を含むより現実的な値，橙は政府統計値

球面調和関数（３）

晩夏・小暑・温風至（あつかぜいたる）

球面調和関数（２）からの続き

あいからわず，球面調和関数に到達していない。Juliaで特殊関数を使うにはどうすればよいか調べる。まあ普通は SpecialFunctions.jl というパッケージを入れればよいと思うのだが，ここには限られたものしかなくて，直交多項式のたぐいは含まれていないので他を探すことになる。ちょっと面倒だ。

この点，Mathematicaはいいのかわるいのか統一的に管理されているので探し回る手間はない。その特殊関数のところをみると，ガンマ関数，ベータ関数，誤差関数，指数積分，直交多項式，ベッセル関数，ルジャンドル関数，双曲線関数，楕円積分，楕円関数，モジュラ形式，ゼータ関数，多重対数関数，マシュー関数，回転楕円体関数，ホイン関数…となっていた。

さて，Juliaである。なるべく多くの特殊関数を１つのパッケージで扱いたいので探したところ，GNU-Scienfitic Library のパッケージ GSL.jl を見つけた。早速試してみると，

using GSL
using Plots
using LinearAlgebra

f1(n,z)=sf_hermite_func(n,z)
x1=0.0:0.01:5.0
y3=[f1.(i,x3) for i=0:4]
g1=plot(x1,y1)

f2(l,z)=sf_legendre_Pl(l,z)
x2=0.0:0.01:1.0
y3=[f2.(i,x2) for i=0:4]
g2=plot(x2,y2)

f3(n,z)=sf_laguerre_n(n,0,z)
x3=0.0:0.1:5
y3=[f3.(i,x3) for i=0:4]
g3=plot(x3,y3)

f4(l,m,z)=sf_legendre_sphPlm(l,m,z)*2*sqrt(pi)
x4=0.0:0.01:1.0
y4=[f4.(i,i-1,x4) for i=1:4]
pushfirst!(y4,f4.(0,0,x4))
g4=plot(x4,y4)

plot(g1,g2,g3,g4,legend=:bottomright)

球面調和関数そのものではないが，球面調和関数の規格化因子を用いたルジャンドル陪関数が用意されているのでほぼそのまま使える。エルミート多項式のところには１次元の量子力学的調和振動子波動関数がそのまま定義されていた。さすがに水素原子波動関数はそうは問屋がおろさない。

図　Juliaでの特殊関数の例

東京オリンピック（３）

東京オリンピック（２）からの続き

東京都議会選で自民・公明が過半数割れして都民ファーストが踏みとどまったことから，東京五輪無観客作戦が進められているような報道がされている。結局は徳仁天皇の意向が微妙に無視しきれなかったからなのかもしれない。

その朝日新聞の説は，収容人員が50%で5000人を超える会場は無観客にするというものだ。ただし，IOC関係者やスポンサー枠など五輪貴族様は例外扱いで，犠牲になる小中学生もカウントにはいれないのだろう。それでも，とりあえず計算してみるとつぎのとおり。

　　　会場　　　　　　　　　　　　　収容定員　日数　　延人数

1 オリンピックスタジアム　　　 68,000 　13 　0

2 東京体育館　　　　　　　　　 7,000 　13 　45,500

3 国立代々木競技場　　　　　　 10,200 　16 　0

4 日本武道館　　　　　　　　　 11,000 　11 　0

5 東京国際フォーラム　　　　　 5,000 　10 　25,000

6 国技館　　　　　　　　　　　 7,300 　15 　54,750

7 馬事公苑　　　　　　　　　　 9,300 　12 　55,800

8 武蔵野の森総合スポーツプラザ 7,200 　11 　39,600

9 東京スタジアム　　　　　　　 48,000 　10 　0

10 武蔵野の森公園　　　　　　　 -

11 有明アリーナ　　　　　　　　 15,000 　16 　0

12 有明体操競技場　　　　　　　 12,000 　14 　0

13 有明アーバンスポーツパーク　 7,000 　  8 　28,000

14 有明テニスの森　　　　　　　 19,900 　  9 　0

15 お台場海浜公園　　　　　　　 5,500 　  5 　13,750

16 潮風公園　　　　　　　　　　 12,000 　15 　0

17 青海アーバンスポーツパーク　 8,400 　  9 　37,800

18 大井ホッケー競技場　　　　　 15,000 　14 　0

19 海の森クロスカントリーコース 16,000 　  1 　0

20 海の森水上競技場　　　　　　 16,000 　14 　0

21 カヌー・スラロームセンター　 7,500 　  6 　22,500

22 夢の島公園アーチェリー場　　 5,600 　  9 　25,200

23 東京アクアティクスセンター　 15,000 　27 　0

24 東京辰巳国際水泳場　　　　　 4,700 　16 　37,600

25 札幌大通公園　　　　　　　　 -

26 幕張メッセAホール・Bホール　 9,000 　20 　90,000

27 釣ヶ崎海岸サーフィンビーチ　 6,000 　  8 　24,000

28 さいたまスーパーアリーナ　　 21,000 　15 　0

29 陸上自衛隊朝霞訓練場　　　　 3,200 　10 　16,000

30 霞ヶ関カンツリークラブ　　　 25,000 　  8 　0

31 江の島ヨットハーバー　　　　 3,000 　11 　16,500

32 伊豆ベロドローム　　　　　　 3,600 　  7 　12,600

33 伊豆 MTBコース　　　　　　　 11,500 　  2 　0

34 富士スピードウェイ　　　　　 22,000 　  3 　0

35 福島あづま球場　　　　　　　 14,300 　  3 　0

36 横浜スタジアム　　　　　　　 35,000 　13 　0

37 札幌ドーム　　　　　　　　　 41,000 　  5 　0

38 宮城スタジアム　　　　　　　 49,000 　  6 　0

39 茨城カシマスタジアム　　　　 40,000 　  8 　0

40 埼玉スタジアム2002　　　　　 64,000 　  8 　0

41 横浜国際総合競技場　　　　　 72,000 　  8 　0

   　　合　　計　　　　　　　          　　　　　　　544,700

合計は，54万人ということになった。これに特別枠が加われば，ざっと100万人というところで，当初の目論見であった300万人の1/3にはなる。1日平均5万人の直接の人出効果に相当する。抽選の手間が省けて便利なのかもしれないが，開会式も閉会式も陸上も水泳もサッカーも野球も「いわゆる無観客」なのか。

球面調和関数（２）

球面調和関数（１）からの続き

ぜんぜん，球面調和関数にはなっていないけれど。JuliaのPlotsの散布図で色をつける方法がなんとなくわかった。各点の色に対応する配列 c を用意して，カラースキームを指定して，おまじないの marker_z = c，をオプションにすればいいようだ。

とりあえず，マーカーをcricleにして半径を関数値にあわせて変化させる。色についてはこれとは独立な関数を用意すればよい。α値なども同様に変えられるかもしれない。

using Plots
using Random
using ColorSchemes

Random.seed!(0)
N=10000

x, y, z, c, t, s = rand(N).-0.5, rand(N).-0.5, rand(N).-0.5, rand(N), rand(N), rand(N)
r = (x.*x .+ y.*y .+ z.*z).^0.5
t = acos.(z./r)
s = atan.(y./x)
c = 100 .* z.^2 .* exp.((-5).*r) .+ 0.5

#println(x[10],y[10],z[10],c[10])

scatter(x, y, z,
marker_z = c,
markershape = :circle,
markersize = c,
markeralpha = 0.5,
markercolor = :jet,
markerstrokewidth = 0.2,
markerstrokealpha = 0.3,
markestrokecolor = :white,markerstrokestyle = :dot,
legend = :right,
camera = (60,60)
)

図 Julia の Plots の散布図の例

国土地理院の空中写真

 国土地理院が，地図・空中写真閲覧サービスを提供している。これはタイムマシンみたいものだ。すべての年代を網羅しているわけではないけれど，今はなくなってしまった風景の一部を見ることができる。

早速，明治から1990年頃までの金沢の実家付近を探してみると，1946.11，1947.11，1948.5，1952.11，1953.4，1962.5，1966.7 ，1973.5，1975.10，1982.10，1983.5，1987.10 などのデータが見つかった。

空中写真については，「本サービスでダウンロード可能な空中写真は，出典の明示等を行っていただければ利用可能です（申請不要）」ということだ。

1962年にはここに住んでいたが（1957-1972），今はなき金沢市立工業高等学校とそのグラウンドが見える。泉野小学校，野田中学校，金沢泉丘高等学校の旧校舎もはっきりとわかる。

写真：1962.5の金沢寺町台（国土地理院の空中写真より）

球面調和関数（１）

tsujimotterさんが日曜化学：量子力学の基本と球面調和関数の可視化というテーマで python のmatplotlib を使ったコードを書いていたところ，東北大学の堀田さんに褒められていた。

なるほど，これをJuliaでやってみようかと思ったのだが，Plots の三次元散布図のカラーオプションの制御方法がわからなくてとりあえず挫折した。

Mathematicaだと簡単に実例が見つかるんだわこれが。

With[{lmax = 2}, Table[If[Abs[m] <= l,
SphericalPlot3D[
Abs[2 SphericalHarmonicY[l, l-m, t, p]]^2, {t, 0, Pi}, {p, 0,
2 Pi}, PlotRange -> {{-0.5, 0.5}, {-0.5, 0.5}, {-0.5, 0.5}},
BoxRatios -> {1, 1, 1}, PlotPoints -> 30, Axes -> False,
Boxed -> False, Mesh -> False, ColorFunctionScaling -> False,
ColorFunction ->
Function[{x, y, z, t, p, r}, Blend[{Blue, Orange}, (Cos[Arg[SphericalHarmonicY[l, l - m, t, p]]] + 1)/
2]]], Null], {l, 0, lmax}, {m, 0, lmax}]] //
GraphicsGrid[#, AspectRatio -> 1] &

図　球面調和関数の例 by Mathematica（筑波大学　武内修さんより引用）

 

失敗知識データベース

 畑村洋太郎さん（1940-）は，失敗学を提起し2002年には特定非営利活動法人の失敗学会を設立している。また，個人では，畑村創造工学研究所の名前でも活動している。

これに関連して，失敗学会には失敗知識データベースのページがある。明治時代から2008年までの様々な分野での失敗事例が，事例詳細，経過，原因，対処，知識化，背景などにわたって丁寧に記録分析されていて，たいへん有益である。その後の更新が止まっているのが残念だ。

福島原発事故が含まれておらず，医療・健康分野がその他にごくわずかな例しか取り上げられていないことも惜しい。せっかくだから，教育やその他の行政分野にわたっても失敗知識を蓄積してほしいところ。2chあがりのネトウヨ担当大臣がデジタル庁にうつつを抜かすくらいならば，失敗庁をつくるほうが日本の将来にとってよほど有益だろう。

日高川入相花王

 WOWOWのシネマ歌舞伎で，玉三郎の日高川入相花王渡し場の段（2006）をやっていた。歌舞伎では，これは人形振り（役者が人形遣いに操られている人形のように振る舞う）で演じられることになっているらしい。15年前なので，清姫の坂東玉三郎も人形遣いの尾上菊之助もともに若い。船頭は，坂東薪車だった市川九團次だ。

17世紀はじめに出雲阿国によってはじまった歌舞伎であるが，現在も演じられている主要な作品の数々は，元禄時代に近松門左衛門，竹田出雲，三好松洛，並木宗輔らの戯作者によって人形浄瑠璃の演目として作られて大流行した作品が歌舞伎に移植されたものだ。歌舞伎の人形振りというのはその過程を反映しているのかもしれない。

日高川入相花王，文楽では何回か見ているが，玉三郎の所作は文楽での人形の動きを非常に丁寧にシミュレートしていて美しい。これがMIKIKOによるPerfumeの振り付けの源流になるのだろうか？

高野山から国道371号線を南に進むと，和歌山県と奈良県の県境にある護摩壇山に至る。日高川はこのあたりから始まって，椿山ダムを経て御坊市の河口まで続いている。河口の北には道成寺がある。安珍を追いかける清姫は日高川に飛び込んで蛇に化け，道成寺まで到達して鐘に隠れた安珍を焼き殺すわけだ。

図：伝土佐光重（土佐派）画『道成寺縁起』より

氷室開き

 NHKの列島ニュースで金沢局から氷室開きがとりあげられていた。冬の間に積もった雪や氷を室に蓄えて，夏にかけて利用する氷室の起源は，奈良県天理市福住町の都祁氷室にあるらしい。金沢でも，江戸時代には湯涌温泉あたりに氷室がいくつか設けられていた。

金沢の氷室からは，加賀藩が江戸の徳川将軍に氷を献上していた。7月1日がその献上する氷の出発の日であり「氷室の日」となっていて，その前日6月30日が氷室開きになっている。もっとも江戸時代は旧暦の六月朔日が氷室の日だったようだ。

いまではどうかわからないが，小学校のころは氷室の日に学校で氷室万頭が配られた。もしかすると半ドンくらいの勢いだったかもしれない。いまみると氷室万頭はけっこうおいしそうなのだけれど，当時はそれほどでもなかったな。

写真：湯涌温泉の氷室の仕込み初め（金沢市観光協会から引用）

P. S. 高木屋の茶色の氷室饅頭が小学校の記憶に最も近い。

［１］金沢ふるさと学習　指導資料（金沢市）

［２］金沢の人はどこの氷室饅頭を食べているのか（高井寧香）