算額（その1164）

算額（その1164）

九八 鴻巣市三ツ木山王 三木神社 明治28年(1895)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円1個，楕円2個，正六角形

正六角形の中に，交差する 2 個の半楕円と，楕円に内接する円を容れる。円の直径が 1 寸のとき，正六角形の一辺の長さはいかほどか。

片方の楕円（赤で描いたもの）の長軸が x 軸上にあるように，算額の図を反時計回りに 60° 回転させて考える。もう一つの楕円（灰色で描いたもの）は考えなくてもよい。

正六角形の中心を原点とする
正六角形の一辺の長さを s
楕円の長半径，短半径 を a, b
楕円と正六角形の辺との接点座標を (x, 0)
内接円の半径を r
とおく。
x = s*cos(π/6)
a = 2x
b = s/2
である。

楕円と内接円に関する算法助術の公式84を解き，正六角形の一辺の長さを求める。

include("julia-source.txt");

using SymPy
@syms s, r, a, b
x = s*cosd(Sym(30))
a = 2x
b = s/2
eq1 = x^2 - (a^2 - b^2)*(b^2 - r^2)/b^2
ans_s = solve(eq1, s)[2];

正六角形の一辺の長さは，内接円の半径の √22/2 倍である（直径の √22/4 倍）。
内接円の直径が 1 寸のとき，正六角形の一辺の長さは √22/4 = 1.1726039399558574 寸である。

ans_s |> println
ans_s(r => 1/2).evalf() |> println

   sqrt(22)*r/2
   1.17260393995586

function draw(r, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   s = √22r/2
   (x, y) = s.*(cosd(30), sind(30))
   (a, b) = (2x, s/2)
   plot([0, -x, -x, 0, x, x, 0], [s, y, -y, -s, -y, y, s], color=:blue, lw=0.5)
   circle(0, 0, r, :green)
   circle(-2x, 0, r, :gray70)
   ellipse(-x, 0, a, b, color=:red)
   ellipse(-x/2, s - y/2, a, b, φ=-60, color=:gray70)
   if more        
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(-x, 0, "楕円:a,b,(-x,0)", :red, :center, delta=-delta/2)
       point(r, 0, "r ", :green, :right, :bottom, delta=delta/2)
       point(x, y, "(x,y) ", :blue, :right, delta=-delta/2)
       point(-2x, 0, "-2x", :gray70, :center, delta=-delta/2)
       point(x, 0, " x", :red, :left, :bottom, delta=delta/2)
       point(0, s, " s", :red, :left, :bottom, delta=delta/2)
       point(0, 0, " 0", :black, :left, :bottom, delta=delta/2)
   end
end;

算額（その1163）

算額（その1163）

六四 加須市不動岡 総願寺 慶応2年(1866)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：楕円，長方形，円，正方形，斜線4本

長方形の中に頂点と辺の中点を結ぶ 4 本の斜線を引き，区画された領域の一つに楕円を容れる。楕円の短径が 1 寸のとき，長方形の短辺の長さはいかほどか。

楕円は円を引き伸ばし（引き縮め）たものである。提示された図の長方形の横幅を縮め，縦横比が 1:1 すなわち正方形になったとき楕円は円になる。

縦方向の比率は変わらないので，円の直径から正方形の一辺を求める問題になる。直角三角形の 3 辺とそれに内接する円の直径の関係式を解けば，正方形の一辺の長さがわかる。

長方形の短辺の長さを a，楕円の短半径を r として，以下の方程式を解く。
注：短径は「短半径の2倍」である。

include("julia-source.txt");

using SymPy
@syms a, r
eq = a + a/2 - sqrt(a^2 + a^2/4) - 2r
res = solve(eq, a)[2]  # 2 of 2
res |> println
res(r => 1/2).evalf() |> println

   r*(sqrt(5) + 3)
   2.61803398874989

長方形の短辺の長さ a は，楕円の短半径 r の (√5 + 3) 倍 である。
楕円の短径が 1 寸のとき，長方形の短辺の長さは 2.61803398874989 寸である。

function draw(倍率, more=false)
   pyplot(size=(倍率*400, 400), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1/2
   a = r*(√5 + 3)
   plot(倍率.*[0, a, a, 0, 0], [0, 0, a, a, 0], color=:blue, lw=0.5)
   plot!(倍率.*[0, a/2, a], [a, 0, a], color=:green, lw=0.5)
   plot!(倍率.*[a, 0, a], [a, a/2, 0], color=:green, lw=0.5)
   ellipse(倍率*r, r, 倍率*r, r, color=:red)
   if more        
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, a, " a", :blue, :left, :bottom, delta=delta/2)
       point(倍率*r, r, "(倍率*r,r)", :red, :center, delta=-delta/2)
       point(倍率*a, 0, "(倍率*a,0)", :blue, :center, delta=-delta/2)
       point(倍率*a/2, 0.9a, @sprintf("倍率 = %g", 倍率), :black, :center, mark=false)
       plot!(ylims=(-8delta, a + 5delta), xlims=(-8delta, 倍率*a + 15delta))
   end
end;

算額（その1162）

算額（その1162）

八六 加須市多聞寺 愛宕神社 明治13年(1880)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円8個，外円，二等辺三角形
#Julia, #SymPy, #算額, #和算

外円の中に二等辺三角形を容れ，甲円 3 個，乙円 4　個を容れる。乙円の直径が 2 寸のとき，甲円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
二等辺三角形の底辺の頂点の座標を (x, y); x = sqrt(R^2 - y^2)
甲円の半径と中心座標を r1, (0, y + r1), (x1, y1)
乙円の半径と中心座標を r2, (x2, y2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

@syms R::positive, x::positive, y::negative, r1::positive, x1::positive,
     y1::positive, r2::positive, x2::positive,
     y2::negative
@syms h, r1, r2, r3, x3, y31, y32, y33
eq1 = dist2(0, R, x, y, 0, y + r1, r1)
eq2 = dist2(0, R, x, y, x1, y1, r1)
eq3 = dist2(0, R, x, y, x2, y2, r2)
eq4 = x1^2 + y1^2 - (R - r1)^2
eq5 = x2^2 + y2^2 - (R - r2)^2
eq6 = (x2 - x1)^2 + (y1 - y2)^2 - (r1 + r2)^2
eq7 = y1/x1 - x/(R - y)
eq8 = x^2 + y^2 - R^2;

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return Float64.(v), r.f_converged
end;

function H(u)
   (R, x, y, r1, x1, y1, x2, y2) = u
   return [
-R^2*r1^2 + R^2*x^2 + 2*R*r1^2*y - 2*R*r1*x^2 - 2*R*x^2*y - r1^2*y^2 + 2*r1*x^2*y + x^2*y^2,  # eq1
-R^2*r1^2 + R^2*x^2 - 2*R^2*x*x1 + R^2*x1^2 + 2*R*r1^2*y - 2*R*x^2*y1 + 2*R*x*x1*y + 2*R*x*x1*y1 - 2*R*x1^2*y - r1^2*x^2 - r1^2*y^2 + x^2*y1^2 - 2*x*x1*y*y1 + x1^2*y^2,  # eq2
-R^2*r2^2 + R^2*x^2 - 2*R^2*x*x2 + R^2*x2^2 + 2*R*r2^2*y - 2*R*x^2*y2 + 2*R*x*x2*y + 2*R*x*x2*y2 - 2*R*x2^2*y - r2^2*x^2 - r2^2*y^2 + x^2*y2^2 - 2*x*x2*y*y2 + x2^2*y^2,  # eq3
x1^2 + y1^2 - (R - r1)^2,  # eq4
x2^2 + y2^2 - (R - r2)^2,  # eq5
-(r1 + r2)^2 + (-x1 + x2)^2 + (y1 - y2)^2,  # eq6
-x/(R - y) + y1/x1,  # eq7
-R^2 + x^2 + y^2,  # eq8
   ]
end;
r2 = 2/2
iniv = BigFloat[19, 10, -16.6, 7, 11.5, 3, 12, -8]
res = nls(H, ini=iniv)

   ([4.266666666666667, 2.0655911179772892, -3.7333333333333334, 1.6, 2.581988897471611, 0.6666666666666666, 2.6334969275741744, -1.9328230761165115], true)

甲円の直径は乙円の直径の 1.6 倍である。
乙円の直径が 2 のとき，甲円の直径は 3.2 である。

その他のパラメータは以下の通りである。

   R = 4.26667;  x = 2.06559;  y = -3.73333;  r1 = 1.6;  x1 = 2.58199;  y1 = 0.666667;  r2 = 1;  x2 = 2.6335;  y2 = -1.93282
   二等辺三角形の斜辺の長さ = 4.26667

function draw(r2, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, x, y, r1, x1, y1, x2, y2) = [19, 9.4, -16.6, 7, 11.5, 3, 12, -8]
   (R, x, y, r1, x1, y1, x2, y2) = res[1]
   @printf("乙円の直径が %g のとき，甲円の直径は %g である。\n", 2r2, 2r1)
   @printf("R = %g;  x = %g;  y = %g;  r1 = %g;  x1 = %g;  y1 = %g;  r2 = %g;  x2 = %g;  y2 = %g\n", R, x, y, r1, x1, y1, r2, x2, y2)
   @printf("二等辺三角形の斜辺の長さ = %g\n", sqrt(x^2 + y^2))
   plot()
   circle(0, 0, R)
   circle2(x1, y1, r1, :blue)
   circle(0, y + r1, r1, :blue)
   circle2(x2, y2, r2, :green)
   plot!([x, 0, -x, x], [y, R, y, y], color=:magenta, lw=0.5)
   θ = atand(y1, x1) + atand(-y2, x2)
   (x22, y22) = transform(x2, y2, deg = 2θ)
   circle2(x22, y22, r2, :green)     
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, y + r1, "甲円:r1,(0,y+r1)", :blue, :center, delta=-delta/2)
       point(x1, y1, "甲円:r1,(x1,y1)", :blue, :center, delta=-delta/2)
       point(x2, y2, "乙円:r2\n(x2,y2)", :green, :center, delta=-delta/2)
       point(0, R, " R", :red, :left, :bottom, delta=delta/2)
       point(R, 0, " R", :red, :left, :bottom, delta=delta/2)
   end
end;

---

甲円の半径 だけを求めるには以下のようにすれば，解析解を求めることができる。
a を二等辺三角形の斜辺の長さとする。

using SymPy
@syms a, r1, r2, R
eq1 = a^2 - 16r1*(R - r1)
eq2 = -4R*r1*(R - r1)/(R - 2r1) + a^2
eq3 = r2 - (R*r1 - r1^2)/R
res = solve([eq1, eq2, eq3], (r1, a, R))[2]

   (8*r2/5, 32*sqrt(15)*r2/15, 64*r2/15)

res[1](r2 => 1).evalf() |> println
res[2](r2 => 1).evalf() |> println
res[3](r2 => 1).evalf() |> println

   1.60000000000000
   8.26236447190916
   4.26666666666667

甲円の半径 r1 は，乙円の半径 r2 の 8/5 = 1.6 倍である。
乙円の直径が 2 寸のとき，甲円の直径は 3.2 寸である。

 

 

算額（その1161）

算額（その1161）

一八 大里郡岡部村岡 稲荷社 文化14年(1817)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円6個，長方形，斜線2本

長方形の中に 2 本の斜線を引き，等円 5 個，白円 1 個と，矢（真ん中の等円の頂部から斜線の交点までの線分）を容れる。矢の長さが 1 寸のとき，白円の直径を求めよ。

斜線の交点を原点とする。
斜線はデタラメに引いたものではなく，長方形の長辺を一辺とする正方形の対角線である（図参照）。

矢を「矢」
白円の半径と中心座標を r1, (0, r1)
等円の半径と中心座標を r2, (0, 矢 - r2)
長方形の長辺，短辺を 2a，a + 矢 + 2r2
とおく。

等円の半径は r2 = 矢 + √2矢 = 矢*(1 + √2) である。
⊿OAB と ⊿OCD は相似で，相似比は 3r2:(矢 + 2r2) である。
白円の半径は r1 = (矢 + √2矢) * 3r2/(矢 + 2r2) である。

矢 = 1 のとき，白円の半径は (矢 + √2矢) * 3r2/(矢 + 2r2) = 3 である。

include("julia-source.txt");

function draw(矢, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 矢 + sqrt(2)*矢
   a = 3r2
   b = a
   y = -2r2 - 矢
   相似比 = 3r2/(矢 + 2r2)
   r1 = r2*相似比
   @printf("r2 = %g;  矢 = %g;  相似比 = %g;  r1 = %g\n", r2, 矢, 相似比, r1)
   plot([a, a, -a, -a, a], [-b, b, b, -b, -b], color=:blue, lw=0.5)
   segment(-a, y, a, y, :blue)
   segment(-a, a, a, -a, :blue)
   segment(a, a, -a, -a, :blue)
   circle(0, y + r2, r2)
   circle2(2r2, y + r2, r2)
   circle2(2r2, y + 3r2, r2)
   circle(0, a - r1, r1, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       dimension_line(0, -矢, 0, 0, "矢")
       point(a, y, " (a,y)", :blue, :left, :vcenter)
       point(a, -b, " (a,-b)", :blue, :left, :vcenter)
       point(a, b, " (a,b)", :blue, :left, :vcenter)
       point(0, y+r2, "等円:r2,(0,y+r2)", :red, :center, delta=-delta/2)
       point(0, b - r1, "白円:r1,(0,R-r1)", :green, :center, delta=-delta/2)
       point(0, 0, "O", :blue, :center, :bottom, delta=delta/2)
       point(a, a, "A", :blue, :center, :bottom, delta=delta/2)
       point(-a, a, "B", :blue, :center, :bottom, delta=delta/2)
       point(y, y, "C", :blue, :center, delta=-delta/2)
       point(-y, y, "D", :blue, :center, delta=-delta/2)
       xlims!(-a - 5delta, a + 10delta)
   end
end;

 

算額（その1160）

算額（その1160）

九八 鴻巣市三ツ木山王 三木神社 明治28年(1895)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円4個，半円1個，斜線2本

半円の中に 2 本の斜線を引き，甲円，乙円，丙円を容れる。丙円の直径が 1 寸のとき，乙円の直径はいかほどか。

半円の半径と中心座標を R, (0, 0); R = 2r1
甲円の半径と中心座標を r1, (0, r1)
乙円の半径と中心座標を r2, (x2, r2); r2 = r1/2
丙円の半径と中心座標を r3, (0, R - r3)
半円の端点から乙円への接線と半円の交点座標を (x0, sqrt(R^2 - x0^2))
とおき，以下の連立方程式を解くが，いくつかのパラメータは簡単な関係式を満たす。
「問」は丙円の直径が既知で，乙円の直径を求めよとなっているが，図を描く手順は以下のようにするのがよい。
まず甲円は外円の 1/2 倍，乙円の 2 倍である。半円の端点から乙円への接線を引き，この接線に接するように r3 を決める。
R = 2r1，r2 = r1/2，x2 = √2r1 なので（eq1, eq3 を解けば求まる），実質求めなければならないのは x0, r1 の 2 つである。

include("julia-source.txt")

using SymPy

@syms R::positive, r1::positive, r2::positive, x2::positive,
     r3::positive, y3::positive, x0::positive
R = 2r1
r2 = r1/2
x2 = √Sym(2)r1
eq1 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2 |> expand
eq3 = 4R^2 - 16R*r2

eq4 = dist2(x0, sqrt(R^2 - x0^2), -R, 0, x2, r2, r2)
eq5 = dist2(x0, sqrt(R^2 - x0^2), -R, 0, 0, R - r3, r3);

二元連立方程式ではあるが，SymPy では解けないので，まず x0 に付いて解き，それをもう一方の式に代入して r1 を求める。

ans_x0 = solve(eq4, x0)[2]
ans_x0 = ans_x0 |> sympy.sqrtdenest |> simplify |> factor
ans_x0 |> println

   42*r1*(-35 + 384*sqrt(2))/12769

eq5 = eq5(x0 => ans_x0) |> simplify;
@syms d
ans_r1 = solve(eq5, r1)[1] |> apart |> factor
ans_r1 |> println

   8*r3*(1 + 2*sqrt(2))/21

ans_x0 = ans_x0(r1 => ans_r1) |> simplify
ans_x0 |> println

   16*r3*(314*sqrt(2) + 1501)/12769

乙円の半径は，甲円の半径の 1/2 である。

ans_r2 = ans_r1/2
ans_r2 |> println

   4*r3*(1 + 2*sqrt(2))/21

乙円の半径 r2 は，丙円の半径 r3 の (4 + 8√2)/21 倍である。
丙円の直径が 1 寸の解き，乙円の直径は 0.7292242142373696 寸である。

「術」，「答」では (√56 + 7)/21 = 0.6896816558832325 寸としているが，誤りであろう。それでは図を描けない（不適切な図になる）。

その他のパラメータは以下の通りである。

   R = 1.45845;  r1 = 0.729224;  r2 = 0.364612;  x0 = 1.21862;  x2 = 1.03128

function draw(r3, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   #(r1, r2, x0, x2) = res[1]
   r1 = 8*r3*(1 + 2*sqrt(2))/21
   x0 = 16*r3*(314*sqrt(2) + 1501)/12769
   R = 2r1
   r2 = r1/2
   x2 = √2r1
   R = 2r1
   @printf("R = %g;  r1 = %g;  r2 = %g;  x0 = %g;  x2 = %g\n", R, r1, r2, x0, x2)
   plot()
   circle(0, 0, R, beginangle=0, endangle=180)
   circle(0, r1, r1, :magenta)
   circle(0, R - r3, r3, :orange)
   circle2(x2, r2, r2, :blue)
   segment(-R, 0, x0, sqrt(R^2 - x0^2), :green)
   segment(R, 0, -x0, sqrt(R^2 - x0^2), :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, r1, "甲円:r1,(0,r1)", :magenta, :center, delta=-delta/2)
       point(x2, r2, "乙円:r2,(x2,r2)", :blue, :center, delta=-delta/2)
       point(0, R - r3, "丙円:r3,(0,R-r3)", :orange, :center, delta=-delta/2)
       point(x0, sqrt(R^2 - x0^2), "(x0,√(R^2-x0^2))", :green, :right, :bottom, delta=delta/2)
       point(R, 0, " R", :red, :left, :bottom, delta=delta/2)
       point(0, R, " R", :red, :left, :bottom, delta=delta/2)
   end
end;

算額（その1159）

算額（その1159）

九八 鴻巣市三ツ木山王 三木神社 明治28年(1895)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円3個，接線，直線上

直線上に 2 個の等円が互いに接して載っている。左側の等円と直線の接点から右側の等円に接線を引く。左側の等円にできる円弧に小円を容れる。等円の直径が 1 寸のとき，小円の直径はいかほどか。

等円の半径と中心座標を r1, (0, r1), (2r1, r1)
小円の半径と中心座標を r2, (x2, y2)
左の等円と小円と接線の接点座標を (x0, y0)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms r1::positive, r2::positive, x2::positive, y2::positive,
     x0::positive, y0::positive
eq1 = x0^2 + (y0 -r1)^2 - r1^2
eq2 = (x0 - x2)^2 + (y0 -y2)^2 - r2^2
eq3 = y0/x0 - 1//2
eq4 = x2^2 + (r1 - y2)^2 - (r1 -r2)^2
eq5 = x0/2r1 - 2r2/r1
res = solve([eq1, eq2, eq3, eq4, eq5], (r2, x2, y2, x0, y0))[1]

   (r1/5, 16*r1/25, 13*r1/25, 4*r1/5, 2*r1/5)

小円の半径 r2 は，等円の半径 r1 の 1/5 である。
等円の直径が 1 寸のとき，小円の直径は 0.2 寸である。

その他のパラメータは以下の通りである。

  r1 = 0.5;  r2 = 0.1;  x2 = 0.32;  y2 = 0.26;  x0 = 0.4;  y0 = 0.2

function draw(r1, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, x2, y2, x0, y0) = r1 .* (1/5, 16/25, 13/25, 4/5, 2/5)
   @printf("等円の半径が %g のとき，小円の半径は %g である。\n", 2r1, 2r2)
   @printf("r1 = %g;  r2 = %g;  x2 = %g;  y2 = %g;  x0 = %g;  y0 = %g\n", r1, r2, x2, y2, x0, y0)
   plot()
   circle(0, r1, r1)
   circle(2r1, r1, r1)
   circle(x2, y2, r2, :blue)
   x01 = 2x2 - x0
   y01 = 2y2 - y0
   abline(0, 0, y01/x01, 0, 1.6r1, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, r1, "等円:r1,(0,r1)", :red, :center, delta=-delta/2)
       point(2r1, r1, "等円:r1,(0,r1)", :red, :center, delta=-delta/2)
       point(x2, y2, "小円:r2,(x2,y2)", :blue, :right, :vcenter, deltax=-8delta)
       point(x0, y0, "(x0,y0)", :red, :left, delta=-delta/2)
       segment(0, 0, 2r1, r1, :gray70)
   end
end;

算額（その1158）

算額（その1158）

十七 大里郡岡部村岡 稲荷社 文化13年(1816)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円10個，外円，直角三角形

外円の中に直角三角形，甲円，乙円，丙円を 1 個ずつ，丁円，戊円，己円を 2 個ずつ容れる。丁円，戊円の直径がそれぞれ 4 寸 3 寸のとき，己円の直径はいかほどか。

算額(その589）を一段階複雑にしたものである。

直角三角形の斜辺は外円の直径である（中心を通る）。
外円の半径と中心座標を R, (0, 0)
線分 AB, CA, CB の長さを AB,  CA, CB とする。
甲円の半径と中心座標を r1, (x1, y1)
乙円の半径と中心座標を r2, (0, r2 - R)
丙円の半径と中心座標を r3, (r3 - R, 0)
丁円の半径と中心座標を r4, (x41, y41), (x42, y42)
戊円の半径と中心座標を r5, (x5, -OB - r5)
己円の半径と中心座標を r6, (-OA/2 - r6, y6)
とおき，以下の連立方程式を解く。
まず eq1, eq2, eq3, eq4 を解いて，r6, R, CA, CB を求める。
eq1, eq2, eq3 は「算法助術」の公式29 である。

include("julia-source.txt")

using SymPy

@syms AB::positive, CA::positive, CB::positive,
     R::positive, r1::positive, r2::positive, r3::positive,
     r4::positive, r5::positive, r6::positive
AB = 2R
eq1 = AB^2 - 16R*r4
eq2 = CA^2 - 16R*r5
eq3 = CB^2 - 16R*r6
eq4 = CA^2 + CB^2 - AB^2
(ans_r6, ans_R, ans_CA, ans_CB) = solve([eq1, eq2, eq3, eq4], (r6, R, CA, CB))[1]

   (r4 - r5, 4*r4, 8*sqrt(r4)*sqrt(r5), 8*sqrt(r4)*sqrt(r4 - r5))

己円の半径 r6 は，丁円の半径 r4 から戊円の半径 r5 を引いたものである。
丁円，戊円の直径がそれぞれ 4 寸 3 寸のとき，己円の直径は 1 寸である。

算額の答えはここまででよいが，図を描くために残りのパラメータをすべて求める。

---

乙円の半径はすぐに求まり，次いで戊円の中心座標 x5 も求まる。

@syms x5
ans_r2 = (ans_R - ans_CB/2)/2
eq5 = x5^2 + (r2 - r5)^2 - (r2 + r5)^2
ans_x5 = solve(eq5, x5)[2]
ans_x5 |> println

   2*sqrt(r2)*sqrt(r5)

同様にして，丙円の半径，己円の中心座標 y6 も求まる。

@syms y6
r3 = (ans_R - ans_CA/2)/2 
eq6 = (r3 - r6)^2 + y6^2 - (r3 + r6)^2
ans_y6 = solve(eq6, y6)[2]
ans_y6 |> println

   2*sqrt(2)*sqrt(r6)*sqrt(-sqrt(r4)*sqrt(r5) + r4)

甲円と丁円のパラメータについては，直角三角形の斜辺を水平にした場合について求め，実際の角度だけ座標を回転させて図を描く。

@syms x4
r1 = ans_R/2
eq7 = x4^2 + r4^2 - (ans_R - r4)^2
ans_x4 = solve(eq7, x4)[2]
ans_x4 |> println

   2*sqrt(2)*r4

function draw(r4, r5, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
  (r6, R, CA, CB) = (r4 - r5, 4*r4, 8*sqrt(r4)*sqrt(r5), 8*sqrt(r4)*sqrt(r4 - r5))
   AB = sqrt(CA^2 + CB^2)
   r2 = (R - CB/2)/2
   x5 = 2*sqrt(r2)*sqrt(r5)
   r3 = (R - CA/2)/2 
   y6 = 2*sqrt(2)*sqrt(r6)*sqrt(-sqrt(r4)*sqrt(r5) + r4)
   r1 = R/2
   (x1, y1) = transform(0, r1, deg=-atand(CB/CA))
   x4 = 2√2r4
   (x41, y41) = transform(x4, r4, deg=-atand(CB/CA))
   (x42, y42) = transform(-x4, r4, deg=-atand(CB/CA))
   @printf("丁円，戊円の直径が %g, %g のとき，己円の直径は %g である。\n", 2r4, 2r5, 2r6)
   @printf("r4 = %g;  r5 = %g;  r6 = %g;  R = %g;  CA = %g;  CB = %g\n", r4, r5, r6, R, CA, CB)
   plot([-CA/2, CA/2, -CA/2, -CA/2], [-CB/2, -CB/2, CB/2, -CB/2], color=:green, lw=0.5)
   circle(0, 0, R, :orange)
   circle(0, r2 - R, r2, :blue, lw=2)
   circle2(x5, -CB/2 - r5, r5, :blue)
   circle(r3 - R, 0, r3, :magenta, lw=2)
   circle22(-CA/2 - r6, y6, r6, :magenta)
   circle(x1, y1, r1, lw=2)
   circle(x41, y41, r4)
   circle(x42, y42, r4)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, R, "R", :orange, :center, :bottom, delta=delta/2)
       point(R, 0, " R", :orange, :left, :bottom, delta=delta/2)
       point(x1, y1, "甲円:r1,(x1,y1)", :red, :center, delta=-delta/2)
       point(x41, y41, "丁円:r4\n(x41,y41)", :red, :center, delta=-delta/2)
       point(x42, y42, "丁円:r4\n(x42,y42)", :red, :center, delta=-delta/2)
       point(0, r2 - R, "乙円:r2\n(0,r2-R)", :blue, :center, :bottom, delta=delta/2)
       point(x5, -CB/2 - r5, "戊円:r5\n(x5,-CB/2-r5)", :blue, :center, :bottom, delta=delta/2)
       point(r3 - R, 0, "   丙円:r3,(r3-R,0)", :magenta, :left, :vcenter)
       point(-CA/2 - r6, y6, "   己円:r6,(-CA/2-r6,y6)", :magenta, :left, :vcenter)
       point(-CA/2, -CB/2, "C ", :green, :right, :vcenter)
       point(CA/2, -CB/2, " A", :green, :left, :vcenter)
       point(-CA/2, CB/2, "B ", :green, :right, :vcenter)
   end
end;

算額（その1157）

算額（その1157）

十七 大里郡岡部村岡 稲荷社 文化13年(1816)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円11個，半円2個，長方形

長方形の中に甲円，乙円，丙円の 11 個，半円 2 個を容れる。甲円の直径が 5 寸 5 分のとき，乙円の直径はいかほどか。

長方形の長辺，短辺を 2r1, r1
甲円の半径と中心座標を r1, (0, 0)
乙円の半径と中心座標を r2, (r1 + r2, 0), (x2, r1 - r2)
丙円の半径と中心座標を r3, (x3, y3)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms r1::positive, r2::positive, x2::positive,
     r3::positive, x3::positive, y3::positive
eq1 = (r1 + r2)^2 + r1^2 - (2r1 - r2)^2
eq2 = x2^2 + (2r1 - r2)^2 - (2r1 + r2)^2
eq3 = x3^2 + (y3 + r1)^2 - (2r1 - r3)^2
eq4 = (r1 + r2 - x3)^2 + y3^2 - (r2 + r3)^2
eq5 = x3^2 + y3^2 - (r1 + r3)^2;
res = solve([eq1, eq2, eq3, eq4, eq5], (r2, x2, r3, x3, y3))[1]

   (r1/3, 2*sqrt(6)*r1/3, 2*r1/11, 12*r1/11, 5*r1/11)

丙円の半径 r3 は，甲円の半径 r1 の 2/11 倍である。
甲円の直径が 5.5 寸のとき，丙円の直径は 1 寸である。

その他のパラメータは以下の通りである。

  r1 = 2.75;  r2 = 0.916667;  x2 = 4.49073;  r3 = 0.5;  x3 = 3;  y3 = 1.25

function draw(r1, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, x2, r3, x3, y3) = (r1/3, 2*sqrt(6)*r1/3, 2*r1/11, 12*r1/11, 5*r1/11)
   @printf("甲円の直径が %g のとき，丙円の直径は %g である。\n", 2r1, 2r3)
   @printf("r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g;  y3 = %g\n", r1, r2, x2, r3, x3, y3)
   plot(r1 .* [2, 2, -2, -2, 2], r1*[-1, 1, 1, -1, -1], color=:green, lw=0.5)
   circle(0, -r1, 2r1, :magenta, beginangle=0, endangle=180)
   circle(0, r1, 2r1, :magenta, beginangle=180, endangle=360)
   circle(0, 0, r1)
   circle2(r1 + r2, 0, r2, :green)
   circle4(x2, r1 - r2, r2, :green)
   circle4(x3, y3, r3, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(2r1, 0, " 2r1", :magenta, :left, :bottom, delta=delta/2)
       point(0, r1, " r1", :magenta, :left, :bottom, delta=delta/2)
       point(0, 0, "甲円:r1,(0,0)", :red, :center, delta=-delta)
       point(r1 + r2, 0, "乙円:r2\n(r1+r2,0)", :green, :center, :bottom, delta=delta)
       point(x2, r1 - r2, "乙円:r2\n(x2,r1-r2)", :green, :center, :bottom, delta=delta)
       point(x3, y3, "丙円:r3,(x3,y3)   ", :blue, :right, :vcenter)
       xlims!(-2r1 - 3delta, 2r1 + 8delta)
   end
end;

2024/07/19

オクラの花です

ハイビスカスと同じ仲間です

算額（その1156）

算額（その1156）

九八 鴻巣市三ツ木山王 三木神社 明治28年(1895)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円7個，外円，半円1個，

外円の中に，天円 2 個，地円 2 個，等円 2個半を入れる。天円の直径が 1 寸のとき，地円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
天円の半径と中心座標を r1, (x1, 4r2 - R)
地円の半径と中心座標を r3, (x3, 2r2 - R)
等円の半径と中心座標を r2, (0, r2 - R), (0, 3r2 - R), (0, 5r2 - R)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms R::positive, r1::positive, x1::positive,
     r2::positive, r3::positive, x3::positive
eq1 = x1^2 + (4r2 - R)^2 - (R - r1)^2
eq2 = x3^2 + (2r2 - R)^2 - (R - r3)^2
eq3 = x1^2 + r2^2 - (r1 + r2)^2
eq4 = x3^2 + r2^2 - (r2 + r3)^2
eq5 = (5r2 - R)^2 + r2^2 - R^2
res = solve([eq1, eq2, eq3, eq4, eq5], (R, x1, r2, r3, x3))[1]

   (39*r1/10, 2*r1, 3*r1/2, 4*r1/3, 2*sqrt(13)*r1/3)

地円の半径 r3 は，天円の半径 r1 の 4/3 倍である。
天円の直径が 1 寸のとき，地円の直径は 4/3 である。

その他のパラメータは以下の通りである。

  r1 = 0.5;  R = 1.95;  x1 = 1;  r2 = 0.75;  r3 = 0.666667;  x3 = 1.20185

function draw(r1, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, x1, r2, r3, x3) = r1 .* (39/10, 2, 3/2, 4/3, 2√13/3)
   @printf("天円の直径が　%g のとき，地円の直径は %g である。\n", 2r1, 2r3)
   @printf("r1 = %g;  R = %g;  x1 = %g;  r2 = %g;  r3 = %g;  x3 = %g\n", r1, R, x1, r2, r3, x3)
   y = 5r2 - R
   plot()
   circle(0, 0, R)
   circle2(x1, 4r2 - R, r1, :blue)
   circle2(x3, 2r2 - R, r3, :green)
   circle(0, r2 - R, r2, :magenta)
   circle(0, 3r2 - R, r2, :magenta)
   circle(0, y, r2, :magenta, beginangle=180, endangle=360)
   x = sqrt(R^2 - y^2)
   segment(-x, y, x, y, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(x1, 4r2 - R, "天円:r1\n(x1,4r2−R)", :blue, :center, delta=-delta/2)
       point(x3, 2r2 - R, "地円:r3,(x3,2r2-R)", :green, :center, delta=-delta/2)
       point(0, r2 - R, "等円,(0,r2-R)", :magenta, :center, delta=-delta/2)
       point(0, 3r2 - R, "等円,(0,3r2-R)", :magenta, :center, delta=-delta/2)
       point(0, 5r2 - R, "等円,(0,5r2-R)", :magenta, :center, delta=-delta/2)
   end
end;

算額（その1155）

算額（その1155）

九八 鴻巣市三ツ木山王 三木神社 明治28年(1895)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円4個，楕円，正方形，菱形

正方形の中に楕円，菱形を 1 個ずつ，等円を 4 個容れる。楕円の長径，短径が 4 寸，3 寸のとき，等円の直径はいかほどか。

正方形の一辺の長さを 2s
楕円の長半径，短半径を a, b
楕円と菱形の一辺との接点座標を (x0, y0)
等円の半径と中心座標を r, (s - r, s - r)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms s::positive, r::positive,
     x0::positive, y0::positive, a::positive, b::positive
eq1 = dist2(s, 0, 0, s, s - r, s - r, r)
eq2 = x0^2/a^2 + y0^2/b^2 - 1
eq3 = -b^2*x0/(a^2*y0^2) + 1
eq4 = y0/(s - x0) - 1
res = solve([eq1, eq2, eq3, eq4], (r, s, x0, y0))[1]

   ((a*sqrt(4*a^2 + 1) - a + sqrt(2)*b*sqrt(sqrt(4*a^2 + 1) - 1) - sqrt(2*a^4 - a^2*sqrt(4*a^2 + 1) + a^2 + sqrt(2)*a*b*sqrt(4*a^2 + 1)*sqrt(sqrt(4*a^2 + 1) - 1) - sqrt(2)*a*b*sqrt(sqrt(4*a^2 + 1) - 1) + b^2*sqrt(4*a^2 + 1) - b^2))/(2*a), (a*sqrt(4*a^2 + 1) - a + sqrt(2)*b*sqrt(sqrt(4*a^2 + 1) - 1))/(2*a), sqrt(4*a^2 + 1)/2 - 1/2, b*sqrt(2*sqrt(4*a^2 + 1) - 2)/(2*a))

@syms t
res[1](sqrt(4a^2 +  1) => t) |> println
res[2](sqrt(4a^2 +  1) => t) |> println
res[3](sqrt(4a^2 +  1) => t) |> println
res[4](sqrt(4a^2 +  1) => t) |> println

   (a*t - a + sqrt(2)*b*sqrt(t - 1) - sqrt(2*a^4 - a^2*t + a^2 + sqrt(2)*a*b*t*sqrt(t - 1) - sqrt(2)*a*b*sqrt(t - 1) + b^2*t - b^2))/(2*a)
   (a*t - a + sqrt(2)*b*sqrt(t - 1))/(2*a)
   t/2 - 1/2
   b*sqrt(2*t - 2)/(2*a)

楕円の長半径 a，短半径 b が与えられたとき，等円の半径は
(a*t - a + sqrt(2)*b*sqrt(t - 1) - sqrt(2*a^4 - a^2*t + a^2 + sqrt(2)*a*b*t*sqrt(t - 1) - sqrt(2)*a*b*sqrt(t - 1) + b^2*t - b^2))/(2*a)
ただし，t = sqrt(4a^2 +  1) と，長々しい式になる。

楕円の長半径，短半径が 4/2, 3/2 のとき，等円の直径は 1.463744764599768 である。

「答」は 「1.461 余り」としている。

a = 4/2
b = 3/2
t = sqrt(4a^2 +  1)
r = (a*t - a + sqrt(2)*b*sqrt(t - 1) - sqrt(2*a^4 - a^2*t + a^2 + sqrt(2)*a*b*t*sqrt(t - 1) - sqrt(2)*a*b*sqrt(t - 1) + b^2*t - b^2))/(2*a)
2r

   1.463744764599768

function draw(a, b, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   t = sqrt(4a^2 +  1)
   r = (a*t - a + sqrt(2)*b*sqrt(t - 1) - sqrt(2*a^4 - a^2*t + a^2 + sqrt(2)*a*b*t*sqrt(t - 1) - sqrt(2)*a*b*sqrt(t - 1) + b^2*t - b^2))/(2*a)
   s = (a*t - a + sqrt(2)*b*sqrt(t - 1))/(2*a)
   x0 = t/2 - 1/2
   y0 = b*sqrt(2*t - 2)/(2*a)
   @printf("等円の直径は %g である。\n", 2r)
   plot([s, s, -s, -s, s], [-s, s, s, -s, -s], color=:magenta, lw=0.5)
   plot!([s, 0, -s, 0, s], [0, s, 0, -s, 0], color=:green, lw=0.5)
   circle4(s - r, s - r, r)
   ellipse(0, 0, a, b, color=:blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(s - r, s - r, "等円:r,(s-r,s-r)", :red, :center, delta=-delta/2)
       point(x0, y0)
   end
end;

Julia で統計解析 2024/7/18 版

Julia で統計解析

2024/7/18 にバージョンアップ

https://r-de-r.github.io/stats/Julia-stats1.html
https://r-de-r.github.io/stats/Julia-stats2.html
https://r-de-r.github.io/stats/Julia-stats3.html
https://r-de-r.github.io/stats/Julia-stats4.html
https://r-de-r.github.io/stats/Julia-stats5.html
https://r-de-r.github.io/stats/Julia-stats6.html
https://r-de-r.github.io/stats/Julia-stats7.html
https://r-de-r.github.io/stats/Julia-stats8.html
https://r-de-r.github.io/stats/Julia-stats9.html

 

第 1 章　Julia の実行環境

1.　必要なファイルをダウンロードする
2.　Julia のインストール
3.　Julia を起動し終了する
4.　作業ディレクトリを変える
5.　Julia の環境設定をする
5.1.　startup.jl ファイルがあるかを確認
5.2.　config がなかった場合
5.3.　startup.jl がある場合
5.4.　変更内容の確認
6.　オンラインヘルプを使う
7.　パッケージを利用する
8.　エディタを使う
8.1.　Atom
8.2.　Jupyter lab
9.　結果を保存する

第 2 章　データの入出力

1.　データをデータフレームへ読み込む
1.1.　Excel のワークシートファイル
1.2.　CSV ファイル
1.2.1.　典型的な CSV ファイル
1.2.2.　一般の CSV ファイル
1.2.3.　CSV.read の引数
1.3. インターネット上の CSV ファイル
1.4.　モニター上に表示された表をデータフレームに読み込む
1.5.　クリップボードにコピーした内容をデータフレームにする
1.6.　文字列行から入力する
1.7.　エンコーディングの違うファイルを読む
1.8.　デリミタで区切られていない固定書式ファイルを読む
2.　データフレームを CSV ファイルに書き出す
2.1.　CSV.write の引数

第 3 章　データフレームの取り扱い

1.　データを使用するための準備
1.1.　既存のデータを使用する
1.2.　自前のデータ
2.　 データフレームの概要
2.1.　データフレームの大きさ
2.2.　データフレームの変数名
2.3.　データフレームの表示
3.　データフレームのコピーは copy()  で
4.　空データフレーム
5.　データフレームの列の参照
6.　データフレームの行の参照
7.　データフレームの列名の変更
8.　列の抽出
9.　指定された列を削除する
10.　行の抽出
11.　行の削除
12.　重複を除き，ユニークな行のみを含むデータフレームを作る
13.　欠損値を含む行か含まない行か
14.　欠損値を含まない行を抽出する
15.　データフレームの列方向連結
16.　データフレームの行方向連結
16.1.　cols = :setequal, cols=:orderequal
16.2.　cols = :intersect
16.3.　cols = :subset
16.4.　cols = :union
17.　データフレームの最終行の次に 1 行追加する
18.　行を繰り返してデータフレームを作る
19.　データフレームに列を挿入する
20.　データフレームの要素から新しいデータフレームを作る
21.　データフレームの各列に関数を施す
22.　ソート（並べ替え）
22.1.　ソートされているかをチェックする
22.2.　ソートの順番を決める
22.3.　並べ替えベクトルを返す
23.　行の包含
23.1.　df1 と df2 のすべての組み合わせを作る crossjoin
23.2.　df1 の行のうち，df2 に含まれない行を抽出する antijoin
23.3.　df1 に，df2 をマージする innerjoin
23.4.　df1 に，df2 をマージする leftjoin
23.5.　df1 に，df2 をマージする outerjoin
23.6.　df2 に，df1 をマージする rightjoin
23.7.　両方に存在する項目のみでマージする semijoin
24.　欠損値のリストワイズ除去
25.　ロングフォーマット（縦長データフレーム）に変換する
26.　ワイドフォーマット（横長データフレーム）に変換する
27.　データフレームをグループ変数に基づいて分割する
27.1. グループ分けされたデータフレームを抽出する
27.2.　グループ化されたデータフレームに関する情報
27.3.　親データフレームを返す
28.　列ごとに関数を適用する
29.　行ごとに関数を適用する
30.　クエリーによる操作例
30.1.　要素に関数を適用する　@map コマンド
30.2.　条件を満たす行を抽出　@filter コマンド
30.3.　データフレームのグループ化　@groupby コマンド
30.4.　データのソート @orderby，@orderby_descenidng，@thenby，@thenby_descending コマンド
30.5.　データフレームのマージ　@groupjoin コマンド
30.6.　データフレームの連結　@join コマンド
30.7.　キーと値のペアを展開　@mapmany コマンド
30.8.　要素を取り出す　@take コマンド
30.9.　要素を捨てる　@drop コマンド
30.10.　重複データを除く　@unique コマンド
30.11.　列の選択　@select コマンド
30.12.　列名の変更　@rename コマンド
30.13.　変数変換　@mutate コマンド
30.14.　欠損値行を除く　@dropna コマンド
30.15.　欠損値を含まないデータフレーム　@dissallowna コマンド
30.16.　欠損値の置き換え　@replacena コマンド
31.　データフレームを二次元配列に変換する

第 4 章　一変量統計

1.　一変量統計
1.1.　一変数の場合
1.1.1.　基礎統計量
1.1.2.　パーセンタイル値
1.1.3.　度数分布
1.2.　複数の変数の場合
1.2.1.　eachcol() を使う
1.2.2.　describe() を使う
1.2.3.　combine(), select()/select!(), transform()/transform!() を使う
1.3.　グループごとの記述統計量
1.3.1.　describe() を使う
1.3.2.　変数ごとに統計量を一覧表示
2.　二変量統計
2.1.　二重クロス集計表
2.2.　相関係数，共分散
2.2.1.　欠損値を含まない場合
2.2.1.1.　それぞれの関数を使う
2.2.1.2.　combine() を使う
2.2.2.　欠損値を含む場合
2.2.2.1.　それぞれの関数を使う
2.2.2.2.　combine() を使う
2.3.　相関係数行列，分散・共分散行列
2.3.1.　欠損値を含まない場合
2.3.2.　欠損値を含む場合
3.　多変量統計
3.1. 多重クロス集計表

第 5 章　離散データの可視化

1.　例示に使用するデータセット
1.1.　カテゴリーデータ
2.　棒グラフ
2.1.　一標本の場合
2.2.　二標本以上の場合
2.2.1.　横に並べる棒グラフ
2.2.2.　積み上げ棒グラフ
2.3.　複数のグラフを行列状にまとめて表示する方法
3.　帯グラフ
4.　モザイクプロット
5.　バルーンプロット

第 6 章　数値データの可視化

1.　ヒストグラム
1.1.　一標本の場合
1.2.　複数標本の場合
2.　ボックスプロット（箱ひげ図）
3.　バイオリンプロット
4.　ドットプロット
5.　カーネル密度推定の描画
6.　Q-Q プロット
7.　散布図
8.　カーネル密度推定
9.　散布図行列
10.　作図レシピ
10.1.　二軸グラフ
10.2.　Andrews plot
10.3.　星座グラフ
10.4.　サンフラワープロット
10.5.　レーダーチャート

第 7 章　検定と推定

1.　検定関数関数の呼び出し方
2.　検定関数関数により得られる結果の利用法
3.　分布の検定
3.1.　観察度数が一様かどうかの検定
3.1.1.　ピアソンの \\chi^2 検定
3.1.2.　対数尤度比検定（G^2 検定）
3.2.　観察度数が理論比に從うかどうかの検定
3.2.1.　ピアソンの \\chi^2 検定
3.2.2.　対数尤度比検定（G^2 検定）
4.　独立性の検定
4.1.　ピアソンの \\chi^2 検定
4.2.　対数尤度比検定（G^2 検定）
5.　パワーダイバージェンス検定
6.　フィッシャーの正確検定
7.　二項検定
8.　t 検定
8.1.　一標本の検定（母平均の検定）
8.2.　等分散の場合の t 検定
8.3.　等分散でない場合 Welch の方法
9.　マン・ホイットニーの U 検定
10.　符号検定
11.　ウィルコクソンの符号付き順位和検定
12.　相関係数の検定
13.　対応のない k 標本（独立 k 標本）
13.1.　一元元配置分散分析
13.2.　一元元配置分散分析（ウェルチの方法）
13.3.　クラスカル・ウォリス検定
14.　コルモゴロフ・スミルノフ検定
14.1.　1 標本の分布の検定
14.2.　2 標本の分布の差の検定
15.　アンダーソン・ダーリング検定
15.1.　1 標本の場合
15.2.　k 標本の場合
16.　ワルド・ウォルフォビッツ連検定
17.　並べ替え検定（無作為検定）

第 8 章　多変量解析

1. 回帰分析
1.1.　線形最小二乗回帰 Linear Least Square
1.1.1.　単回帰分析
1.1.2.　重回帰分析
1.2.　リッジ回帰 Ridge Regression
1.2.1.　説明変数が 1 個の場合
1.2.2.　説明変数が 2 個以上の場合
1.3.　GLM パッケージによる回帰分析
1.3.1.　重回帰分析 Ordinary Least Squares Regression
1.3.2.　プロビット回帰 Probit Regression
1.3.3.　ロジット回帰 Logit Regression
1.4.　多項式回帰
1.4.1.　重回帰プログラムを用いる
1.4.2.　多項式回帰 Polynomials パッケージ
1.5.　指数回帰
1.5.1.　説明変数が 1 個の場合
1.5.2.　説明変数が 2 個以上の場合
1.6.　累乗回帰
1.6.1.　説明変数が 1 個の場合
1.6.2.　説明変数が 2 個以上の場合
1.7.　非線形回帰
2.　判別分析
2.1.　二群判別分析
2.2.　多重判別分析
3.　主成分分析
4.　因子分析
5.　古典的多次元尺度解析
6.　クラスター分析
6.1.　K-means 法による非階層的クラスター分析
6.2. 階層的クラスター分析
7.　カテゴリー変数の取り扱い方
7.1.　重回帰分析の場合
7.2.　判別分析の場合

第 9 章　関数

 

 

算額（その1154）

算額（その1154）

一〇一 大宮市高鼻町 氷川神社 明治31年(1898)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円7個，楕円1個

外円の中に楕円 1 個，甲円 2 個を容れ，隙間に乙円 2 個，丙円 2　個を容れる。乙円の直径が 1 寸のとき，丙円の直径はいかほどか。

注：これだけでは条件不足である。「問」には書かれていないが，乙円は楕円の曲率円である。そのように解釈すれば条件は足りる。

楕円の長半径，短半径と中心座標を a, b, (0, 0); a = r1 + r3, b = a - 2r3
外円の半径と中心座標を R, (0, 0); R = a
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (R - r2, 0)
丙円の半径と中心座標を r3, (0, R - r3)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms a, b, R, r1, r2, r3

a = r1 + r3
b = a - 2r3
R = a
eq1 = (R - r2)^2 + (R - r1)^2 - (r1 + r2)^2
eq2 = r2 - b^2/a
res= solve([eq1, eq2], (r3, r1))[3]  # 3 of 3

   (r2*(sqrt(2) + 2)/2, r2*(4 + 3*sqrt(2))/2)

丙円の半径 r3 は，乙円の半径 r2 の √2/2 + 1 倍である。
乙円の直径が 1 寸のとき，丙円の直径は 1.7071067811865475 寸である。

「術」は簡約化すれば，「乙円の径の √2/2 + 1 倍が丙円の径である」となり，同じである。

乙径 = 1
天 = √8 + 3
(天 - √天)*乙径/2

   1.7071067811865475

(√Sym(8) + 3 - sqrt(√Sym(8) + 3))/2 |> sympy.sqrtdenest |> simplify |> println

   sqrt(2)/2 + 1

その他のパラメータは以下の通りである。

  r2 = 0.5;  r3 = 0.853553;  r1 = 2.06066;  a = 2.91421;  b = 1.20711

function draw(r2, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   #@printf("正方形の一辺の長さが %g，大円の直径が %g のとき，小円の直径は %g である。\n", a, 2r1, 2r2)
   (r3, r1) = (r2*(sqrt(2) + 2)/2, r2*(4 + 3*sqrt(2))/2)
   a = r1 + r3
   b = a - 2r3
   @printf("乙円の直径が %g のとき，丙円の直径は %g である。\n", 2r2, 2r3)
   @printf("r2 = %g;  r3 = %g;  r1 = %g;  a = %g;  b = %g\n", r2, r3, r1, a, b)
   plot()
   circle(0, 0, a, :orange)
   circle22(0, a - r1, r1)
   circle22(0, a - r3, r3, :blue)
   circle2(a - r2, 0, r2, :green)
   ellipse(0, 0, a, b, color=:magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, a - r1, "甲円:r1,(0,a-r1)", :red, :center, delta=-delta/2)
       point(a - r2, 0, "乙円:r2\n(a-r2,0)", :green, :center, :bottom, delta=delta/2)
       point(0, a - r3, "丙円:r3,(0,a-r3)", :blue, :center, delta=-delta/2)
       point(0, b, "b", :magenta, :center, delta=-delta/2)
       point(0, a - 2r3, "a-2r3", :blue, :center, :bottom, delta=delta/2)
       point(0, a, "a", :orange, :center, :bottom, delta=delta/2)
       point(a, 0, " a", :orange, :left, :bottom, delta=delta/2)
   end
end;

算額（その1153）

算額（その1153）

　七九 春日部市 香取神社 明治9年(1876)
　九九 春日部市小渕 観音院 明治30年(1897)
一〇三 春日部市 東福寺 明治40年(1907)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：正三角形，正方形，三角法

面積が 65996035.98832495593 坪の正三角形の土地の一辺の長さはいかほどか。

注：10間四方=1坪

正三角形の一辺の長さが a，高さが h のとき，正三角形の面積 S は，S = (a/2) * (√3a/2) = a^2 * √3/4 である。

√3/4 = 0.4330127018922193 であるが，この近似値 0.433 が「三角法」と呼ばれていた。
これは，辺の長さが同じ正三角形の面積(S)と正方形の面積(a^2)の比である。

a = sqrt(S/(√3/4)) ≒ sqrt(S/0.433) 

演算精度が必要なので，長精度演算を行う。

正三角形の面積 = big"65996035.98832495593"

   6.59960359883249559299999999999999999999999999999999999999999999999999999999999e+07

正方形の一辺の長さ = 正三角形の一辺の長さ = sqrt(正三角形の面積/0.433)

   12345.6789000000000569783430598987438553058946472557173681712465830204994722506

三角法の近似値として 0.433 を使うと，正三角形の面積が 65996035.98832495593 坪なので，正三角形の一辺の長さは sqrt(65996035.98832495593/0.433) = 12345.6789 = 123456.789 間である。
この「きれいな結果」を出すために，あのとんでもない長い数字列の面積が与えられたのだ。

sqrt(65996035.98832495593/0.433) |> println
sqrt(big"65996035.98832495593"/big"0.433") |> println

   12345.678899999999
   12345.6788999999999999999999999999999999999999999999999999999999999999999999999

0.433 という近似値ではなく，正確な値 √3/4 を 使うとどうなるか，逆算してみる。
正三角形の一辺の長さが 12345.6789 になる，面積を求める。

一辺 = big"12345.6789"
面積 = 一辺/2 * √big"3" * 一辺/2

   6.599797195723032842573567192825322468159425832771741027437062249938522329618199e+07

三角法 = √3/4
面積 = big"65997971.957230328"
一辺 = sqrt(面積/三角法)

   12345.67890000000031782004072795357343003813226134859822153538170864063729519095

面積/一辺^2

   0.43301270189221929829415103085921145975589752197265625

正確な三角法を使うときには，「正三角形の面積は 65997971.957230328 坪」といえばよい。そうすれば答えは 123456.789 間になる。

算額（その1152）

算額（その1152）

　九九 春日部市小渕 観音院 明治30年(1897)
一〇三 春日部市 東福寺 明治40年(1907)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：正五角形，正方形，五角法

面積が 262155154.5032769612 坪の正五角形の土地の一辺の長さはいかほどか。

以下の図は GeoGebra で，一辺を与えた正多角形をかき，その面積を求めたものである。算額の「答」を得るために，「問」の数値を変えている。

注：10間四方=1坪

正五角形の一辺の長さを a，一辺の両端と正五角形の中心を結ぶ二等辺三角形の高さを h とする。
h = a/2/tand(Sym(72)/2) で，正五角形の面積は，5 * a*h/2 = 5 * (a * (a/2)/tand(72/2))/2 である。

include("julia-source.txt")

using SymPy

@syms S::positive, a::positive

S = 5 * a*h/2
S |> println

   5*a^2/(4*sqrt(5 - 2*sqrt(5)))

a^2 にかかる係数 5/(4*sqrt(5 - 2*sqrt(5))) = 1.7204774005889671 を「五角法」と呼ぶ。

5/(4*sqrt(5 - 2*sqrt(5)))

   1.7204774005889671

慣例的に，五角法の近似値として 1.72 が使われていた。その場合，正五角形の面積は S = 1.72*a^2 となる。a について解けば a = sqrt(S/1.72) である。

与えられた条件は，正五角形の面積が 262155154.5032769612 なので，正五角形の一辺の長さは sqrt(262155154.5032769612/1.72) = 12345.6789 = 123456.789 間である。
この「きれいな結果」を出すために，あのとんでもない長い数字列の面積が与えられたのだ。

sqrt(262155154.5032769612/1.72) |> println
sqrt(big"262155154.5032769612"/big"1.72") |> println

   12345.6789
   12345.67890000000000000000000000000000000000000000000000000000000000000000000004

しかし，近似値でない五角法を使えば，「きれいな数値は出てこない」

sqrt(262155154.5032769612/1.7204774005889671) |> println
sqrt(big"262155154.5032769612"/big"1.7204774005889671") |> println

   12343.96593263111
   12343.96593263111017155386297511751038450862787310648018338764402293502103615974

最初の図に示したが，正確な五角法を使うときには，「正五角形の面積は 262227917.89 坪」といえばよい。そうすれば答えは 123456.789 間になる。

sqrt(big"262227917.89"/big"1.7204774005889671")

   12345.67890000004107475558293899828913755690970317874535819580788216955810780953