算額（その1574）

算額（その1574）

三重県鳥羽市 観世音 寛政12年(1800)
深川英俊，トニー・ロスマン：聖なる数学：算額，p. 97，森北出版株式会社，2010年４月22日.

キーワード：円1個，二等辺三角形2個
#Julia, #SymPy, #算額, #和算

外円の中に大小の正三角形を容れる。外円の直径が与えられたとき，小正三角形の一辺の長さはいかほどか。

外円の半径と中心座標を R, (0, 0)
小正三角形の一辺の長さを q
とおき，以下の方程式を解く。

include("julia-source.txt")
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R, q
eq = (q/2)^2 + (-R/2 - √Sym(3)*q/2)^2 - R^2
res = solve(eq, q)[1]
res |> println

    R*(-sqrt(3) + sqrt(15))/4

小正三角形の一辺の長さ q は，外円の半径 R の (√15 - √3)/4 倍である。
外円の直径が 100 のとき，小正三角形の一辺の長さは 26.76165673298175 である。

100/2*(√15 - √3)/4

    26.76165673298175

function draw(R, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    q = R*(√15 - √3)/4
    @printf("外円の直径が = %g のとき，小正三角形の一辺の長さは %g である。\n", 2R, q)
    plot()
    circle(0, 0, R)
    polygon(0, 0, R, 3, color=:blue)
    polygon(0, -R/2 - q/√3, q/√3, 3, 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(0, R, "R", :red, :center, :bottom, delta=delta/2)
        point(0, -R/2, "-R/2", :blue, :center, :bottom, delta=delta/2)
        point(√3R/2, -R/2, "(√3R/2,-R/2)  ", :blue, :right, :bottom, delta=delta/2)
        point(q/2, -R/2 - √3q/2, "(q/2,-R/2-√3q/2)", :blue, :left, delta=-delta/2)
    end  
end;

draw(100/2, true)

 

算額（その1573）

算額（その1573）

岡山県瀬戸内市長船町土師宮森 片山日子神社 明治6年(1873)
http://www.wasan.jp/okayama/katayamahiko.html
深川英俊，トニー・ロスマン：聖なる数学：算額，p. 96，森北出版株式会社，2010年４月22日.

キーワード：円1個，二等辺三角形
#Julia, #SymPy, #算額, #和算

二等辺三角形に円を容れる。底辺と斜辺が与えられたとき，円の直径はいかほどか。

底辺と斜辺を a, b
円の半径と中心座標を r, (0, r)
とおき，以下の連立方程式を解く。

include("julia-source.txt")
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a, b, r, h
h = sqrt(b^2 - (a/2)^2)
eq1 = ((a + 2b)*r/2) - (h*a/2)
res = solve(eq1, r)[1]
res |> println

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

円の半径は a*sqrt(4b^2 - a^2)/(2a + 4b) である。
底辺が 12，斜辺が 10 のとき，内接円の直径は 6 である。

2*res(a => 12, b => 10) |> println

    6

以下のほうが簡単。

@syms a, b, r, h
h = sqrt(b^2 - (a/2)^2)
eq2 = (a/2)/b - r/(h - r)
res2 = solve(eq2, r)[1]
res2 |> println

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

function draw(a, b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r = a*sqrt(4b^2 - a^2)/(2a + 4b)
    @printf("a = %g;  b = %g;  r = %g\n", a, b, r)
    plot([a/2, 0, -a/2, a/2], [0, sqrt(b^2 - (a/2)^2), 0, 0], color=:blue, lw=0.5)
    circle(0, r, r)
    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, sqrt(b^2 - a^2/4), "sqrt(b^2-a^2/4)", :blue, :center, :bottom, delta=delta/2)
        point(a/2, 0, "a/2", :blue, :left, :bottom, delta=delta/2)
        point(0, r, " r", :red, :left, :vcenter)
    end  
end;

draw(12, 10, true)

大島の「純手打」さぬきうどん 大島うどん

高松市太田上町 大島の「純手打」さぬきうどん 大島うどん

Julia の歪度と尖度

1.　Julia の歪度と尖度は古典的定義に基づく

1.1.　skewness, kurtosis は StatsBase に含まれる。

# Skewness
# This is Type 1 definition according to Joanes and Gill (1998)
"""
    skewness(v, [wv::AbstractWeights], m=mean(v))

Compute the standardized skewness of a real-valued array `v`, optionally
specifying a weighting vector `wv` and a center `m`.
"""

# (excessive) Kurtosis
# This is Type 1 definition according to Joanes and Gill (1998)
"""
    kurtosis(v, [wv::AbstractWeights], m=mean(v))

Compute the excess kurtosis of a real-valued array `v`, optionally
specifying a weighting vector `wv` and a center `m`.
"""

1.2.　R は 3 通りの定義に基づく関数値を返す

Joanes and Gill (1998) [^1]discuss three methods for estimating kurtosis and skewness.

[1]: D. N. Joanes and C. A. Gill (1998), Comparing measures of sample skewness and kurtosis. The Statistician, 47, 183–189.

1.2.1.　Kurtosis

Type 1: This is the typical definition used in many older textbooks.

g2 = m4/m2^2 - 3

Type 2: Used in SAS and SPSS.

G2 = ((n + 1)*g2 + 6)*(n − 1)/((n − 2)*(n − 3))

Type 3: Used in MINITAB and BMDP.

b2 = m4/s^4 − 3 = (g2 + 3)*(1 − 1/n)^2 − 3

Only G2 (corresponding to type = 2) is unbiased under normality.

デフォルトは type=3 である。

1.2.2.　Skewness

Type 1: This is the typical definition used in many older textbooks.

g1 = m3/m2^(3/2)

Type 2: Used in SAS and SPSS.

G1 = g1*sqrt(n*(n - 1))/(n - 2)

Type 3: Used in MINITAB and BMDP.

b1 = m3/s^3 = g1*((n - 1)/n)^(3/2)

All three skewness measures are unbiased under normality.

デフォルトは type=3 である。

1.3.　Julia でも type を選べる関数を定義する

1.3.1.　尖度

using StatsBase

kurt1(x) = kurt(x, type=1)
kurt2(x) = kurt(x, type=2)
kurt3(x) = kurt(x)
function kurt(x; type=3)
    doc = """
    kurt(x; type=n)
    type: an integer between 1 and 3 selecting one of the algorithms for computing kurtosis detailed below.

    Type 1: g2. This is the typical definition used in many older textbooks.
    Type 2: G2. Used in SAS and SPSS.
    Type 3: b2. default. Used in MINITAB and BMDP.
    """
    type in 1:3 || error(doc)
    n = length(x)
    g2 = kurtosis(x)
    if type == 2
        return  ((n + 1)*g2 + 6)*(n − 1)/((n − 2)*(n − 3))
    elseif type== 3
        return (g2 + 3)*(1 − 1/n)^2 − 3
    else
        return g2
    end
end;           

x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 1];

kurt(x), kurt3(x)

    (-1.6335639958376689, -1.6335639958376689)

kurt(x, type=1), kurt1(x)

    (-1.3130419701699616, -1.3130419701699616)

kurt(x, type=2), kurt2(x)

    (-1.356984911550468, -1.356984911550468)

kurt(x, type=3), kurt3(x)

    (-1.6335639958376689, -1.6335639958376689)

1.3.2.　歪度

using StatsBase

skew1(x) = skew(x, type=1)
skew2(x) = skew(x, type=2)
skew3(x) = skew(x)
function skew(x; type=3)
    doc = """
    skew(x; type=n)
    type: an integer between 1 and 3 selecting one of the algorithms for computing skewness detailed below.

    Type 1: g1. This is the typical definition used in many older textbooks.
    Type 2: G1. Used in SAS and SPSS.
    Type 3: b1. default. Used in MINITAB and BMDP.
    """
    type in 1:3 || error(doc)
    n = length(x)
    g1 = skewness(x)
    if type == 2
        return g1*sqrt(n*(n - 1))/(n - 2)
    elseif type== 3
        return g1*((n - 1)/n)^(3/2)
    else
        return g1
    end
end;           

x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 1]
skew(x), skew3(x)

    (0.10905343709973031, 0.10905343709973031)

skew(x, type=1), skew1(x)

    (0.12772490663150174, 0.12772490663150174)

skew(x, type=2), skew2(x)

    (0.15146310708295876, 0.15146310708295876)

skew(x, type=3), skew3(x)

    (0.10905343709973031, 0.10905343709973031)

 

算額（その1572）

算額（その1572）

四十 岩手県一関市牧沢 牧沢八幡神社 明治5年(1872)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

一関市博物館 第23回 令和6年度　和算に挑戦！
https://www.city.ichinoseki.iwate.jp/index.cfm/7,177295,c,html/177295/20241125-090413.pdf

キーワード：円，長方形，折り紙
#Julia, #SymPy, #算額, #和算

長辺，短辺が 4，3 の長方形の紙を折ってできる直角三角形に内接する円の直径はいかほどか。

山村の図は相変わらず不明瞭・不正確で（そもそも何を描いているのかさえわからない），算額の穴埋めも失敗している。「一関博物館」は流石である。

長辺，短辺の長さを a, b
折り線と長辺の交点座標を (c, 0)
元の頂点が対角線上に重なる座標を (x0, y0)
円の半径と中心座標を r, (x, y)
とおき，以下の連立方程式を解く。

円の半径だけを求めるならば eq1, eq2 だけの連立方程式を解けばよい。

include("julia-source.txt")
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a, b, c, r, diag, x0, y0, x, y
eq1 = ((sqrt(a^2 + b^2) - b)^2 + c^2) - (a - c)^2
eq2 = c + sqrt(a^2 + b^2) - b - (a - c) - 2r
eq3 = y0/(a - x0) - b/a
eq4 = (x0 - c)^2 + y0^2 - c^2;

res = solve([eq1, eq2, eq3, eq4], (c, r, x0, y0))[1]

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

対角線の長さを d = sqrt(a^2 + b^2) とすれば，
折り線と長方形の頂点の交点の x 座標は c = b*(d - b)/a
円の半径は (d - a - b)/2 + c である。

以下は図を描くために円の中心座標を求める連立方程式であるが，a, b を記号のままで解くと長い式になるので，eq5, eq6 の a, b に特定の値を与えて解くほうが現実的である。

eq5 = dist2(a, 0, 0, b, x, y, r)(r => res[2], x0 => res[3], y0 => res[4])(a => 4, b => 3)
eq6 = dist2(c, 0, x0, y0, x, y, r)(c => res[1], r => res[2], x0 => res[3], y0 => res[4])(a => 4, b => 3);

solve([eq5, eq6], (x, y))[3]  # 3 of 4

    (5/2, 1/2)

function draw(a, b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    d = sqrt(a^2 + b^2)
    c = b*(d - b)/a
    r = (d - a - b)/2 + c
    (x0, y0) = (a*b/d, -b^2/d + b)
    (x, y) = (5/2, 1/2)
    @printf("a = %g;  b = %g;  d = %g;  c = %g;  r = %g\n", a, b, d, c, r)
    plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:green, lw=0.5, seriestype=:shape, fillcolor=:lemonchiffon)
    plot!([c, a, x0, c], [0, 0, y0, 0], color=:gray, seriestype=:shape, fillcolor=:orange)
    plot!(a .- [c, a, x0, c], b .- [0, 0, y0, 0], color=:gray, seriestype=:shape, fillcolor=:orange)
    plot!([0, c, x0, 0], [b, 0, y0, b], color=:red, seriestype=:shape, fillcolor=:white)
    plot!(a .- [0, c, x0, 0], b .- [b, 0, y0, b], color=:red, seriestype=:shape, fillcolor=:white)
    circlef(x, y, r)
    circlef(a - x, b - y, r)
    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, b, " b", :blue, :left, :bottom, delta=delta/2) 
        point(a, b, "(a,b)", :blue, :right, :bottom, delta=delta/2) 
        point(a - c, b, "(a-c,b)", :blue, :center, :bottom, delta=delta/2) 
        point(c, 0, "c ", :blue, :right, :bottom, delta=delta/2) 
        point(x0, y0, "(x0,y0)", :blue, :left, :bottom, delta=delta/2)
        point(x, y, "(x,y)", :blue, :left, :bottom, delta=delta/2)
    end  
end;

draw(4, 3, true)