円の長さ

正方形の長さからの続き

ある図形の大きさの指標となる長さを，図形内に一様分布する2点の距離の期待値として定義することで，都道府県の形や大きさを，面積や周長だけでなく"長さ"で特徴づけるという話をしている。

練習として，正方形内のランダムな2点の距離の期待値$\ d \ $を計算できることを確認した。次にトライするのが円であるが，ネットで検索しても生成AIにきいてもあまり適切な解答が得られない。一番近いのが，Yahoo知恵袋の「半径１の円内の任意の２点間の距離の期待値は？」だ。これも結局解析的な答えがでなくて，数値計算で $\ d=0.9054\ $という値を出している。

$d = \int_0^{2\pi} (\frac{1}{2\pi}) d\theta_1 \int_0^{2\pi} ( \frac{1}{2\pi}) d\theta_2 \int_0^1 (2 r_1) dr_1 \int_0^1 (2 r_2) dr_2 \sqrt{r_1^2+r_2^2-2r_1 r_2 \cos(\theta_1-\theta_2) }$

ただし，各積分の()内がそれぞれの変数に対応する確率密度関数，$p(\theta_1)$，$ p(\theta_2)$，$q(r_1)$， $q(r_2)\ $であり，それぞれの変数で積分すると1になるように規格化されている。

これを計算するためには，与えられた4変数の確率分布関数から変数変換によって，積分可能な形に持ち込む必要があるが，なかなか難渋する。しかたがないので，とりあえずJuliaとMathematicaで数値計算してみる。

a=zeros(Float64,1000001,2)

function ju(a,n)

  k = 0

  for i in 1:n

      x = 2*rand()-1

      y = 2*rand()-1

      if x^2+y^2 < 1.0

        k = k + 1

        a[k,1] = x

        a[k,2] = y

      end

  end 

  return k

end

function su(a,n)

  m = div(n,4)*3

  sum = 0

  for i = 1:m

    for j = i:m

      sum = sum + sqrt((a[i,1]-a[j,1])^2+(a[i,2]-a[j,2])^2)

    end

  end

  return sum/binomial(m,2)

end

n=300000

@time su(a,n)

235821
24.223278 seconds (7.18 k allocations: 500.875 KiB, 0.07% compilation time)
0.9055523

生成AIの2つが答えた，$\dfrac{4}{\pi}=1.27324$はたぶん誤っていたということだろう。

解析的に計算できないかと思うのだが，角度積分が完全楕円積分の形になるので，これをさらに積分するのはちょっと難しそうだった。角度積分を後回しにしてもさらに面倒か。Mathematicaに投げてみたが，忍耐可能時間内には答えが出なかった。

都道府県の長さ

ACTIVE GALATTICさんが，次のように書いていた。

中核都市の金沢と被害が激しい輪島は東京と水戸くらい離れていると石川県の広さを説明する投稿を先日見かけたが、石川県の広さは34~35位とランキング下位なので、むしろ東京や大阪が都道府県としては異様に狭いと認識すべきなのだろう

面積と距離が微妙に絡まっていて，ややモニョる。ある図形が与えられたときに，その面積ではなくて，長さを図形の大きさの指標として与えることはできるだろうか。地図上の領域だと東西南北の端点を考えて，南北間距離，東西間距離を抽出できるが，これでは2つの数字の組になる。あるいは図形の周長も考えられるが，フラクタル的な図形だと面積が等しくても周長はいくらでも長くできる。

そこで考えたのが次のような指標だ。図形を格子状の正方形に分割する。これらの正方形の全ての対を考え，各対の正方形の中心間の距離の平均値を求める。分割する正方形のサイズを小さくした極限として，一つの長さが定まるのではないか。（方法１）

検索してみると，もっとスマートな方法があった。図形内の点をランダムに選ぶ。その際に縦横方向にそれぞれ座標値が一様分布であるとした確率分布を仮定する。こうして選んだ2点間の距離の期待値を求めればよいというものだ。（方法２）

まず，方法１のプログラムをJuliaで作成して計算してみた。

a=zeros(Int64,1000001,2)

function ju(a,m,n)

  k = 0

  for i in 1:m

    for j in 1:n

      k = k + 1

      a[k,1] = i

      a[k,2] = j

    end

  end

end

function su(a,m,n)

  mn = m*n

  sum = 0

  for i = 1:mn

    for j = i:mn

      sum = sum + sqrt((a[i,1]-a[j,1])^2+(a[i,2]-a[j,2])^2)

    end

  end

  return sum/binomial(mn,2)

end

m=800

n=800

ju(a,m,n)

@time su(a,m,n)/m

n=   25    time= 0.016720　0.521829

n=   50　time= 0.018188　0.521513

n= 100　time= 0.065998　0.521433

n= 200　time= 0.788281　0.521412

n= 400　time= 12.1924  　0.521407

n= 800　time= 198.615  　0.521405

nは正方形の1辺の格子点の数，timeは計算時間(秒)であり，最後が計算結果の平均長である。

方法２による理論値が解析的に得られており，( 2+√2+5 log(√2 + 1) )/15 = 0.521405 なので，nを増やすとともに理論値に収束し，n=800では理論値に一致している。

m = 50; n = 50;

a = Flatten[Table[{i, j}, {i, 1, m}, {j, 1, n}], 1];

Timing[Sum[

     Sqrt[(a[[i, 1]] - a[[j, 1]])^2 + (a[[i, 2]] - a[[j, 2]])^2], 

　　　{i, 1, n*n}, {j, i, n*n}]/Binomial[n*n, 2]/n // N]

Mathematicaでは，若干コードが簡単になるけれど，時間は n=50 で 12.4801秒と，600倍以上かかってしまった。

（付）このコードはそのまま乱数による計算に置き換えることもできる。Juliaでは，a[k,1] = rand()，a[k,2] = rand() として最後の結果の/nを除く。Mathematicaでは，Tableの中身を{i,j}から{RandomReal[1], RandomReal[1]}に置き換えて最後の結果の/nを除く。

［１］地域内距離（腰塚武志）

JuliaのPyCallとSymPy

Maxima-Jupyterからの続き

Mathematicaの代替となるフリーの数式処理言語として，Maximaが候補にはなったが，やはり何となく使いにくい。Jupyter上のmaximaで終端記号のコロンを忘れるだけで毎回大量のエラーメッセージが出るのにも閉口する。

とすれば，Juliaで数式処理ということになる。標準的な手法としては，Pythonの数式処理モジュールであるSymPyを使えばよいが，先日からPython 呼び出しのPyCallとSymPyがエラーになっている。実害がなかったので放置していたけれど，再度見直してみた。

これまでは，googleで問題解決を図っていたが，必要な情報にたどり着くまではいくつものサイトを渡り歩く必要があった。そこに登場したChatGPTはかなりよい仕事をしてくれるので助かる。問題は，PyCallにあったので，エラーコードを与えて尋ねてみた。

Ｑ：User

Juliaで，PyCallとSymPy がうまく動きません。

例えば，

（中略）

Error building `PyCall`: 

どうしたらいいでしょう。

Ａ：ChatGPT

問題は、おそらく libpython が見つからないために発生しているようです。以下の手順で問題を解決してみてください。

Julia REPLを終了して、再度起動してください。

次に、次のコマンドを実行して、PyCallを再構築し、環境変数を設定してください。

julia

Copy code

ENV["PYTHON"] = "/Users/koshi/.julia/conda/3/bin/python"

import Pkg

Pkg.build("PyCall")

もしまだ問題が解決しない場合、次の手順を試してください。

（中略）

もしまだ問題が解決しない場合、PYCALL_DEBUG_BUILD 環境変数を設定して、デバッグ情報を取得してみてください。

（中略）

これで、エラーが発生している理由に関する追加情報が得られるはずです。その情報を使用して、問題を特定し、解決策を見つけることができます。

SymPyの問題が解決しない場合、まずPyCallが正常に機能していることを確認してください。その後、SymPyの再構築を試みてください。

これらの手順のいずれかが、PyCallとSymPyの問題を解決するはずです。それでも問題が解決しない場合は、エラーメッセージとともに質問してください。

で最初のPYTHON環境変数を /opt/homebrew/bin/python3.11 に置き換えて（だったかな？） なんとかクリアすることができた。

最小交換硬貨枚数（１）

大学入学共通テスト（４）からの続き

2025年度の大学入学共通テストに新規参入する「情報Ｉ」の試作問題が，日経の朝刊に掲載されていた。

退職前には，監督に当たるのはこれ以上体力的精神的な限界をはるかにこえてもう絶対無理と思っていた大学入学センター試験だ。これが，大学入学共通テストと看板を替えて，教科「情報」を新しい試験科目として追加採用してしまった。情報教育関係の大学や高校の教員の皆さんはお喜びのようであるが，なんだかなぁの案件である。

60分で全問必答の大問が4問出題されている。掲載されていた第3問がプログラミングの問題だったので，さっそくチャレンジしてみた。その問題のテーマは最小交換硬貨枚数だ。最近は，PayPayで支払うことが増えてきたが，財布の硬貨数の制御は老人の認知症防止のために最適である。おつりが増えすぎないように，少し硬貨を加えて切りのよいおつりにするあれね。

例えば，46円の支払いだと，10円×４ + 5円×１ + 1円×１ でもよいし，50円×１ + 1円×１ - 5円×１ もある。ここで，マイナスはお店から戻ってくる硬貨を表わしている。前者の交換硬貨枚数は6枚であり，後者では3枚となり，後者の方が交換硬貨枚数は少ない。

問題の誘導部分を読む前に，さっそくプログラムを作ってみたが，肝腎の多めに払っておつりが返ってくるところに不備がありまくりだった。

function pay(m,y)

    n = 0

    c=[500,100,50,10,5,1]

    d=[y]

    for i in 1:6

        while y-c[i] >= 0

            y = y - c[i]

            n = n + 1

            push!(d,c[i])

        end

    end

    push!(d,-(m+n))

    println(d)

    return m+n

end

function change(y)

    m = 0

    n = 0

    c=[1,5,10,50,100,500]

    d=[y]

    pay(0,y)

    for i in 1:5

        while mod(y,c[i+1])!=0

            y = y + c[i]

            n = n + 1

            push!(d,c[i])

        end

        m1 = m + n

        if m1!=m

            push!(d,-m1)

            print(d)

            pay(m1,y)

            n=0

            deleteat!(d,length(d))

        end

        m = m1

    end

end

change(46)

[46, 10, 10, 10, 10, 5, 1, -6]
[46, 1, 1, 1, 1, -4][50, 50, -5]
[46, 1, 1, 1, 1, 50, -5][100, 100, -6]
[46, 1, 1, 1, 1, 50, 100, 100, 100, 100, -9][500, 500, -10]

そこで，1000円以下なら数もしれているということで，方針を変えて総当たりで確認するアルゴリズムに変更した。全くスマートではないのである。この中から最小値を探すのは目の子でできる。出力される配列 [ ] の中身は，交換される硬貨枚数が並んでいて，マイナスがついたものは店から客へのバック分を表わしている。最後の出力値が交換硬貨枚数だ。

function foop(y)
    c=[1,5,10,50,100,500]
    k=[0,0,0,0,0,0]
    for k[1] in -4:4
    for k[2] in -1:1
        for k[3] in -4:4
        for k[4] in -1:1
            for k[5] in -4:4
            for k[6] in 0:2
                z = 0
                n = 0
                for j in 1:6
                    n = n + abs(k[j])
                    z = z + c[j]*k[j]
                end
                if z==y
                    println(z," : ",k," : ",n)
                end
            end
            end
        end
        end
    end
    end
end

foop(46)

46 : [-4, 0, 0, -1, -4, 1] : 10
46 : [-4, 0, 0, -1, 1, 0] : 6
46 : [-4, 0, 0, 1, 0, 0] : 5
46 : [1, -1, 0, -1, -4, 1] : 8
46 : [1, -1, 0, -1, 1, 0] : 4
46 : [1, -1, 0, 1, 0, 0] : 3
46 : [1, 1, -1, -1, -4, 1] : 9
46 : [1, 1, -1, -1, 1, 0] : 5
46 : [1, 1, -1, 1, 0, 0] : 4
46 : [1, 1, 4, 0, 0, 0] : 6

ここまでくるのに2時間くらいかかったので，とてもじゃないけれど自分の場合60分で大問4問も解けそうにはない。ループ変数を配列にするなんていうのは初めての経験だった。


写真：日本の通常硬貨（造幣局から引用）

Julia 1.7.0

11月30日にJulia 1.7.0がリリースされた。普段は，しばらく様子を見てからインストールするのだけれど，macOS Apple Siliconにネイティブ対応したようなので，早速試してみる。

メジャーバージョンアップの際には，いつも困ってしまうのだが，今回もすんなりとはいかなかった。とりあえず手順を書いてみるが，何度もやり直しているので，ヒステリシスがあるかもしれない。

　(1) julia 1.7.0 packege ダウンロード＆インストール
　(2) .zshrcのパス変更 % cat .zshrc
export PATH=/Applications/Julia-1.7.app/Contents/Resources/julia/bin
　(3) ln -s /Applications/Julia-1.7.app/Contents/Resources/julia/bin/julia julia
　(4) julia -> using Pkg -> Pkg.add(“IJulia”) -> Pkg.build(“IJulia”)
　(5) brew install jupyter
　(6) jupyter kernetspec list -> jupyter kernetspec uninstall julia1.6.0
　(7) jupyter notebook -> OK
　(8) pip3 uninstall jupyterlab -> pip3 install jupyterlab
　(9) jupyter lab build -> OK

とりあえず，jupyter notebook も jupyter lab も動いたが，パッケージの整合性などからか，これまでのコードがすべてそのまま動くというわけにはいかなかった。ボチボチ様子見よう。

P. S. とりあえず，Plots 周りの不整合はなんとかなった。たいした処理をしていなくても，ときどきカーネルが死んでしまうのが問題だ。

ルジャンドル変換（１）

 ルジャンドル変換は，積分を使って表現するとわかりやすいという説と通常の説明の対応について考えてみたい。この過程でtikzにおける塗りつぶし方法を練習した。

図：ルジャンドル変換のための説明図

原点を通る単調増加関数上の点C $(x,p)$があって，図のように矩形の領域をとる。矩形領域内で関数の下の部分の面積を$f(x)$，上の部分の面積を$f^*(p)$とすると，$f(x)+f^*(p)=xp$が成り立つ。$f^*(p)$が$f(x)$のルジャンドル変換になる。同様に，$f(x)$は$f^*(p)$のルジャンドル変換である。$f^*(p)$のルジャンドル変換は$f^{**}(x)$とも書けるから，$f^{**}(x)=f(x)$となって，ルジャンドル変換を2回繰り返すと元の関数に戻る。

この単調増加関数は，$x$の関数とすると$f'(x)$であり，$p$の関数と見れば${f^{*}}'(p)$となる。そこで，$f^*(p)$を求めるには，$p=f'(x)$を$x$について解いて，$x=\varphi(p)$を求めてから，$f^*(p) = x p - f(x) = \varphi(p) \cdot p - f(\varphi(p))$として求めることができる。

いいかえれば，$f(x) = x p - f^*(p) |_{p=f'(x)}$ として，元の関数$f(x)$が，傾き$p$と切片$ - f^*(p)$でも表現されるということになる。

／begin{tikzpicture}
\tikzstyle{every node}=[font = \Large];
\filldraw (0,0) circle(1pt) node[below left]{O};
\draw[->] (-2,0) -- (8,0) node[right]{$x$};
\draw[->] (0,-2) -- (0,8) node[above]{$p$};
\draw[step=1.0, dotted] (-2,-2) grid (8,8);
\draw [dotted, pattern=north west lines, pattern color=blue] (0,0) -- (1,1.3) -- (2,2) -- (3,2.4) -- (4,3)-- (5,4) -- (6,5.4) -- (6,0);
\draw [dotted, pattern=north west lines, pattern color=red] (0,0) -- (1,1.3) -- (2,2) -- (3,2.4) -- (4,3) -- (5,4) -- (6,5.4) -- (0,5.4);
\draw[domain=0:3, thick] plot(\x,{-0.2*(\x-3.5)^2+2.45});
\draw[domain=3:7, thick] plot(\x,{0.2*(\x-2)^2+2.2});
\filldraw (6,0) circle(1pt) node[below]{$x$};
\filldraw (3,2.4) circle(1pt);
\filldraw (6,5.4) circle(1pt) node[right]{C};
\node[below right] at (6.5,5.4) {$f^{'}(x)$};
\node[above left] at (6,5.4) {$f^{*'}(p)$};
\filldraw (0,5.4) circle(1pt) node[left] {$p$};
\draw[blue] (4,1) node{$f(x)$};
\draw[red] (2,4) node{$f^*(p)$};
\end{tikzpicture}

日本人のノーベル賞

日本における最近の科学や経済の衰退傾向から，将来，日本人の自然科学系のノーベル賞受賞者が出なくなるのではといわれることがある。一方，最近の中国や韓国は日本より科学技術，経済分野で先んじてはいるけれど，まだノーベル賞受賞者をどんどん輩出するようにはなっていない。

ある研究や発明がなされたときと，それが評価されて普及するようになるまでには時差があることが，その一因と考えられる。そこで，日本人（日本出身）の自然科学系のノーベル賞受賞者25名の，研究・発明年と受賞年をグラフに書いてみた。研究発表の時点から受賞までには平均で約26年かかっているようだ。

図： 日本人の受賞年（横軸）と研究発表年（縦軸）
（赤は時差が25.9年，オレンジは25.9±5年を表す）

湯川秀樹が最初に受賞してからの50年間では5名だけだったのが，2000年以後の20年間余で20名ということになる。1970年代から1990年代にかけて，日本の経済や科学を取り巻く環境が良かった時代を反映しているのかもしれない。

なお，上記のグラフを書くためのJuliaのコードは以下の通りであり，Gadfly.jl でグラフをオーバレイする方法がわかった。

using Gadfly
using Compose
using DataFrames
X = [1949,1965,1973,1981,1987,2000,2001,2002,2002,2008,2008,2008,2008,2010,2010,2012,2014,2014,2014,2015,2015,2016,2018,2019,2021]
Y = [1935,1947,1957,1952,1976,1976,1987,1985,1987,1962,1973,1973,1960,1977,1979,2006,1985,1985,1992,1997,1996,1992,1992,1985,1967]
Labels = ["湯川","朝永","江崎","福井","利根川","白川","野依","田中","小柴","下村","小林","益川","南部","根岸","鈴木","山中","赤碕","天野","中村","大村","梶田","大隅","本庶","吉野","真鍋"]

p1=layer(x=X, y=Y, label=Labels, Geom.point, Geom.label, Theme(major_label_font="CMU Serif",minor_label_font="CMU Serif",major_label_font_size=12pt,minor_label_font_size=12pt))
p2=layer(x->x-25.9, 1949,2035, color=[colorant"red"])
plot(p1,p2)
p3=layer(x->x-20.9, 1949,2035, color=[colorant"orange"])
plot(p1,p2)
p4=layer(x->x-30.9, 1949,2035, color=[colorant"orange"]) 

plot(p1,p2,p3,p4)

MCMCへの道（３）

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

次は，MCMC法#3マルコフ連鎖である。昼食のメニューが３種類に限られていて，なおかつ当日のメニュー選択の確率が，前日のメニューだけによって定まるという例があげられていた。はい，この例は非常によくわかる。

問題は，Juliaの配列と行列の扱いだ。MahthematicaはすべてListと考えればよいのでわかりやすかったが，Juliaでは型として，ArrayもMatrixもある。見様見真似でとりあえずコードを書いたのだけれど，まだ十分理解できていない。このため，3x3行列の配列であるbをベクトルvにかけて得られる配列の転置ができなかったので，手動で転置している。

#using LinearAlgebra
using Plots

x=[zeros(3) for i in 1:11]
y=[zeros(11) for i in 1:3]
#
a=[0.2 0.1 0.3 ; 0.2 0.6 0.5 ; 0.6 0.3 0.2]
b=[zeros(3,3) for i in 1:11]
b[1]=[1 0 0 ; 0 1 0 ; 0 0 1]
v=[0.3,0.2,0.5]
#
for i in 2:11
  b[i]=a*b[i-1]
end

for j in 1:11
  x[j]=b[j]*v
end
#
for i in 1:11
  for j in 1:3
    y[j][i]=x[i][j]
  end
end
#
plot(y,legend=:bottom, xlim=(0,10),ylim=(0,0.6))

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 による棄却法の例

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くらいかかってしまう。

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

コラッツ予想

 コラッツ予想に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までのシーケンス数の散布図

球面調和関数（４）

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

多くの疑問はよく探してみると答えになるウェブサイトがみつかるものだ。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)

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

球面調和関数（３）

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

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

あいからわず，球面調和関数に到達していない。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での特殊関数の例

球面調和関数（２）

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

ぜんぜん，球面調和関数にはなっていないけれど。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 の散布図の例

コロナワクチン接種率

 日本政府（各省庁等）の統計データの管理技術にはなかなか残念な綻びが散見される。コロナワクチン接種率もその一つだ。政権の見栄えのために厚生労働省から引き剥がして内閣府に設けた「新型コロナワクチンについて」のページには，医療従事者等と高齢者等の総ワクチン接種回数のデータが，日別と都道府県別で上がっている。しかし，データは必ずしも定時に定常的にアップロードされるわけではなく，目詰まりしながら間欠的に提示されている。都道府県別データに至っては，5月23日の時点で5月14日や16日のデータが最新という有様だ。

WHOの新型コロナウイルス感染者数と死亡者数でも多くの国では毎日更新データが提供されているが，スペインなど１週間で２回ほど更新される程度の国もある。なお，日本は，このデータについては毎日更新ができている。

この直近のワクチン接種回数の都道府県別データを，各都道府県の人口で除したものを医療従事者と高齢者を２つの要素として散布図にしてみた。juliaのグラフィックスのgadfly.jlを用いると簡単にできることがわかったので，それを使ってみた。

- - - - - - - - - - - - - - - -
using Gadfly
using Compose
using DataFrames
X = [4.0,4.2,5.2,4.1,4.8,4.8,4.4,4.0,4.1,4.2,3.5,3.2,3.6,3.0,4.3,5.0,5.0,6.3,3.6,4.3,4.3,3.7,3.9,4.6,4.0,4.0,3.4,3.4,4.6,5.4,5.6,5.3,4.8,4.6,5.8,6.8,4.8,5.6,6.7,4.0,7.1,5.9,6.0,5.9,5.3,5.7,4.6]
Y = [0.6,1.3,1.3,0.7,2.1,2.1,1.3,0.7,0.7,1.0,0.4,0.4,0.5,0.4,1.1,0.4,1.8,1.4,1.3,1.2,1.0,0.5,0.8,0.7,1.4,0.5,0.5,0.5,1.0,3.5,2.3,1.5,0.4,0.8,2.5,1.6,1.8,0.6,3.2,0.4,1.7,0.9,1.0,1.2,1.8,1.4,0.9]
Labels = ["北海道","青森","岩手","宮城","秋田","山形","福島","茨城","栃木","群馬","埼玉","千葉","東京","神奈川","新潟","富山","石川","福井","山梨","長野","岐阜","静岡","愛知","三重","滋賀","京都","大阪","兵庫","奈良","和歌山","鳥取","島根","岡山","広島","山口","徳島","香川","愛媛","高知","福岡","佐賀","長崎","熊本","大分","宮崎","鹿児島","沖縄"]

plot(x=X, y=Y, label=Labels, Geom.point, Geom.label, Theme(major_label_font="CMU Serif",minor_label_font="CMU Serif",major_label_font_size=12pt,minor_label_font_size=12pt))
#plot(x=X, y=Y, Guide.annotation(compose(context(), text(X, Y, Labels))))
- - - - - - - - - - - - - - - -

図　都道府県別コロナワクチン接種率（x: 医療従事者 %，y: 高齢者 %）

ラベルの重複処理ができているといえばできているのだけれど，欠落しているといえばそうなわけで，これを解決すると（#でコメントアウトしたcomposeによる方法）見にくくなってしまうので一長一短なのだった。

エビデンスベースの政策立案（EBPM）が鳴り物入りで喧伝されているが，公文書を堂々と改竄・隠蔽して憚らない文化の中で，恣意的なエビデンスのチェリーピッキングが横行することを考えれば，地味でも着実に正規化された標準データを蓄積して情報公開を進めることの方が，つまらないデジタル庁を設けるよりもまず最優先の課題だと思われる。

ラベル付き散布図

Rをjupyterでからの続き 

ラベル付き散布図と格闘している。

(1) Excel が一番簡単そうで，調べるとすぐに事例が見つかった。しかし，これはWindows版の場合であり，残念ながらMac版（Microsoft Office 2019 Home & Student）にはその機能がなかった。

(2) Rは統計グラフがお手の物のはずなので多分あるだろう。そこでRをjupyter上に持ってきた。これも簡単に事例が見つかった。見本通りにやってみると，

#今回使うライブラリ

library(maptools)

#データの準備

CC <- c("CAN","DNK","FIN","FRA","DEU","ITA","JPN","KOR","NLD","POL","SWE","GBR","USA")

PAR <- c(13.7,15.7,16.6,16.6,20.2,20.1,22.1,10.3,14.9,13.5,17.6,15.6,12.8)

SE <- c(3.8,7.4,8.6,11.5,8.5,12.2,9.3,1.9,5.4,8.9,9.4,6.2,5.5)

sample <- data.frame(PAR, SE)

rownames(sample) <- CC

以下のerrorが出て（これは例題に対しては実害はなかったが）日本語も化けてしまう。

Checking rgeos availability: FALSE
Note: when rgeos is not available, polygon geometry computations in maptools
depend on gpclib, which has a restricted licence. It is disabled by default;
to enable gpclib, type gpclibPermit()

ただ，日本語文字化けを避けるには（エラーは残るものの），次のようにすればよかった。

par(family= "HiraKakuProN-W3")

plot(sample$PAR, sample$SE, xlab="65歳以上人口比率(%)", ylab="対GDP比高齢者向け社会支出(%)")

pointLabel(x=sample$PAR, y=sample$SE, labels=rownames(sample))

図 R on jupyter のラベル付き散布図の例

(3) juliaではラベル付き散布図の例が探せなかったが，labelled scatter をキーワードとするといとも簡単に解が見つかってしまった。こちらは日本語もすんなり通ってくれる。

角運動量代数（２）

角運動量代数（１）からの続き

Mathematicaの関数には，かなり昔からClebshGordan係数が含まれていた。さすがに素粒子物理を専攻していたウルフラムだと思った。Mathematicaには次の関数が実装されている。

ClebschGordan[{j1,m1},{j2,m2},{j,m}]　|jm〉から|j1m1〉|j2m2〉への分解に対する，クレプシュ・ゴルダン係数を与える。
ThreeJSymbol[{j1,m1},{j2,m2},{j3,m3}]　ウィグナーの3-jシンボルの値を与える（m1+m2+m3=0）。
SixJSymbol[{j1,j2,j3},{j4,j5,j6}]　ラカーの6-jシンボルの値を与える。

Juliaではどうかというと，WingerSymbols.jlというパッケージがあり，以下の関数が使える。
wigner3j(j1,j2,j3,m1,m2,m3=-m1-m2)，wigner6j(j1,j2,j3,j4,j5,j6)，clebschgordan(j1,m2,j2,m2,j3,m3=m1+m2)， racahV(j1,m2,j2,m2,j3,m3=-m1-m2)，racahW(j1,j2,J,j3,J12,J23)

いずれにせよ9-j symbolは実装されていないのだった。一方，PythonのSymPyにあるWignerSymbolsには，Wigner 3-j，6-j，9-j，Clebsch-Gordan係数，Racah係数，Gaunt係数が含まれている。Gaunt係数という名前は知らなかったが，換算行列要素の形では頻出する係数だ。それより，JuliaのWignerFamilies.jlというパッケージで，3-j Symbolの配列が計算されてグラフ化されているのが印象的だった。 

using WignerFamilies
# wigner3j for all j fixing j2=100, j3=60, m2=70, m3=-55, m1=-m2-m3
w3j = wigner3j_f(100, 60, 70, -55)
js = collect(eachindex(w3j))
plot(js, w3j.symbols)

図　f(j) = w3j(j,100,60,-15,70,55) （JuliaPackageから引用）

逆順と平方（３）

逆順と平方（２）からの続き

9桁までの範囲で考えるとしても，３億個のデータのうち条件を満足するのが高々7000個なので，ほとんど無駄な計算をしていることになる。そこで，プログラムは汚くとも，もう少し計算を簡略化できないかと考えた結果，12桁までの計算が数分程度で可能になった。2〜3桁分をまとめてループをまわしているが，この範囲でよいという証明はしたわけではなく，直感に頼っているのだった。

function reverse(n,m)
  nun=zeros(Int,m)
  nre=0
  for k in 1:m
    nun[k] = n % 10
    nre = nre + nun[k]*10^(m-k)
    n = div(n - nun[k],10)
  end
  return nre
end

function counter(n)
  a = [0 1 2 3 10 11 12 13 20 21 22 23 30 31 100 101 102 103 110 111 112 113 120 121 122 123 130 131 200 201 202 203 210 211 212 213 220 221 222 223 230 231 300 301 302 303 310 311 312 313]
  b = [50 4 14]
  c = [10^3 10 10^2]
  cnt=1
  kk = c[n%3+1]
  k1max = b[n%3+1]
  k2max = 1 + div(n+10,14)*49
  k3max = 1 + div(n+10,17)*49
  k4max = 1 + div(n+10,20)*49
  for k1 in 1:k1max
    for k2 in 1:k2max
      for k3 in 1:k3max
        for k4 in 1:k4max
          k=a[k1]+a[k2]*kk+a[k3]*kk*10^3+a[k4]*kk*10^6
          if(k!=0)
            m=Int(ceil(log10(k)+10^(-9)))
            kj=BigInt(k)^2
            mj=Int(ceil(log10(kj)+10^(-9)))
            kr=reverse(k,m)
            krj=BigInt(kr)^2
            if(reverse(kj,mj)==krj)
# println(k,":",kr,":",kr^2,":",reverse(k^2,mj))
              cnt=cnt+1
            end
          end
        end
      end
    end
  end
  return cnt
end

for i in 1:12
  @time println(i," : ",counter(i))
end

1 : 4
0.000491 seconds (176 allocations: 6.000 KiB)
2 : 13
0.000551 seconds (765 allocations: 22.258 KiB)
3 : 37
0.002469 seconds (3.68 k allocations: 96.281 KiB)
4 : 100
0.005823 seconds (18.99 k allocations: 466.195 KiB)
5 : 253
0.021409 seconds (81.69 k allocations: 1.886 MiB)
6 : 615
0.091992 seconds (346.74 k allocations: 7.780 MiB)
7 : 1434
0.441439 seconds (1.61 M allocations: 35.324 MiB, 11.21% gc time)
8 : 3244
1.616126 seconds (6.39 M allocations: 138.140 MiB, 11.31% gc time)
9 : 7116
6.780644 seconds (25.58 M allocations: 546.322 MiB, 13.04% gc time)
10 : 15200
28.310428 seconds (113.32 M allocations: 2.337 GiB, 13.25% gc time)
11 : 23284
108.389549 seconds (437.13 M allocations: 8.934 GiB, 13.71% gc time)
12 : 31368
417.966111 seconds (1.71 G allocations: 34.604 GiB, 14.30% gc time)

逆順と平方（２）

逆順と平方（１）からの続き

中嶋慧さんのpythonによる逆順関数がスマートだったので Juliaでも試してみた。数値を文字列に変換して，逆順にしたものをもういちど数値に戻すというものだ。プログラムとしてはスッキリしたけれど，実行時間が倍になってしまった。

function rev(n)
  b = string(n)
  c = b[end:-1:1]
  d = parse(Int,c)
  return d
end

function counter(m,j)
  cnt = 1
  kmax = Int(ceil(10^(m-1)*sqrt(10)))
  for k in 10^(m-1)+1:kmax
    kr = rev(k)
# kj = BigInt(k)^j
# krj = BigInt(kr)^j
    kj = k*k
    krj = kr*kr
    if(rev(kj)==krj)
# if(2000<=k<=2299)
# println(k,":",kr,":",kr^j,":",rev(k^j))
# end
      cnt = cnt + 1
    end
  end
  return cnt
end

function list(l,j)
  sum=1
  for i in 1:l
    @time sum = sum + counter(i,j)
    println(i," **",j," = ",sum)
  end
end

list(9,2)
0.000010 seconds (36 allocations: 1.781 KiB)
1 **2 = 4
0.000015 seconds (264 allocations: 13.062 KiB)
2 **2 = 13
0.000138 seconds (2.60 k allocations: 128.844 KiB)
3 **2 = 37
0.001734 seconds (25.96 k allocations: 1.254 MiB)
4 **2 = 100
0.016847 seconds (259.48 k allocations: 12.538 MiB)
5 **2 = 253
0.224788 seconds (2.59 M allocations: 125.376 MiB, 12.69% gc time)
6 **2 = 615
2.393325 seconds (25.95 M allocations: 1.224 GiB, 12.89% gc time)
7 **2 = 1434
26.113517 seconds (259.47 M allocations: 12.244 GiB, 11.55% gc time)
8 **2 = 3244
284.531202 seconds (2.59 G allocations: 128.882 GiB, 11.40% gc time)
9 **2 = 7116

逆順と平方（１）

 中嶋慧さんが，次のようなプログラムについての話題を提供していた。ある数nの逆順に並べた数をrev(n)とする。例えば，rev(123)は321である。このとき，rev(n^2)=rev(n)^2を満たすm桁の数nは，何個あるかという問題だ。pythonでは時間がかかるので，juliaにするとよいのではないかという提案だったが，残念ながら，工夫に乏しい10^8個の計算はやたら時間を食うだけだった（追伸：条件が rev(rev(k)^2) = k^2 であると誤解していたので，その後訂正したところ，中嶋さんの結果と一致した）。

function reverse(n,m)
  nun = zeros(Int,m)
  nre = 0
  for k in 1:m
    nun[k] = n % 10
    nre = nre + nun[k]*10^(m-k)
    n = div(n - nun[k],10)
  end
  return nre
end

function counter(m)
  cnt = 1
  kmax = Int64(ceil(10^(m-1)*sqrt(10)))

  for k in 10^(m-1)+1:kmax

    kr = reverse(k,m)

    m2 = Int(ceil(log10(k*k)))

    if(reverse(k*k,m2)==kr*kr)

# if(2000<=k<=2299)
# println(k,":",kr,":",kr^2,":",reverse(kr^2,m2),":",k^2)
# end
      cnt = cnt+1
    end
  end
  return cnt
end

function list(l)

  sum = 0
  for i in 1:l
    @time sum = sum + counter(i)

    println(i," : ",sum)
  end

end

list(9)
0.000003 seconds (6 allocations: 576 bytes)
1 : 4
0.000009 seconds (44 allocations: 4.469 KiB)
2 : 13
0.000097 seconds (434 allocations: 50.859 KiB)
3 : 37
0.000901 seconds (4.33 k allocations: 540.750 KiB)
4 : 100
0.011719 seconds (43.25 k allocations: 5.939 MiB)
5 : 253
0.145435 seconds (432.46 k allocations: 62.688 MiB, 12.05% gc time)
6 : 615
1.100731 seconds (4.32 M allocations: 692.869 MiB, 17.40% gc time)
7 : 1434
11.497505 seconds (43.25 M allocations: 7.088 GiB, 16.11% gc time)
8 : 3244
132.445766 seconds (432.46 M allocations: 77.329 GiB, 16.48% gc time)
9 : 7116