数学検定 過去問題 1級（大学程度・一般） 問題2. cos(６i) －i sin(６i) を求める

数学検定　過去問題
https://www.su-gaku.net/suken/support/past_questions/

1級（大学程度・一般）

〔１級〕１次：計算技能検定

問題2. 次の値を求めなさい。ただし，sin(z) ，cos(z) はそれぞれ複素数 z に関する正弦関数，余弦関数で，i は虚数単位を表します。
 cos(６i) －i sin(６i)

単純にやると結果が数値で出てしまう。

using SymPy
simplify(cos(6im)-complex(0,1)*sin(6im)) # 403.4287934927351 - 0.0im
対数をとってみると
log(cos(6im)-im*sin(6im)) # 6.0 - 0.0im　==> 6
なので
ℯ^6 # 403.4287934927351
であることがわかる。
しかし，これは結果論だな。

Julia の SymPy では，変数の型を指定できる。a, b を複素数として，

@syms a::complex b::complex
simplify(simplify(cos(a) - b * sin(a)).subs(b, im)) # 1.0*exp(-I*a)
となる。
a = 6im
なので
-im * a　# 6 + 0im
よって，1.0ℯ^-ia == ℯ^6

ただし，これも結果がわかっているときに試行錯誤的に解を求めていることに違いはない。
たとえば，
simplify(cos(a) - b * sin(a)).subs(a, 6im)
では
-201.713157370279*I*b + 201.715636122456
になってしまう。

数学検定 過去問題 1級（大学程度・一般） 問題1. ４次方程式

数学検定　過去問題
https://www.su-gaku.net/suken/support/past_questions/

1級（大学程度・一般）

〔１級〕１次：計算技能検定

問題1. ４次方程式

 x^4＋x^3＋k*x^2＋８*x＋６７２０＝０ （kは定数）
の４個の解のうち，２個の解について積を求めたところ８０となりました。このとき，kの値を求めなさい。

using SymPy

@syms a b c d k f g x h3 h2 h1 h0 hab
f = x^4 + x^3 +k*x^2 + 8x + 6720
# 解を -a, -b, -c, -d とすると
g = (x+a)*(x+b)*(x+c)*(x+d)
# 展開して x の次数でまとめる
collect(expand(g), x)
# x^4 + (a+b+c+d)x^3 + (ab+ac+ad+bc+bd+cd)x^2 + (abc+abd+acd+bcd)x + abcd
# f の係数と比較して 5 本の等式を立てる
h3 = Eq(a+b+c+d, 1)
h2 = Eq(a*b+a*c+a*d+b*c+b*d+c*d, k)
h1 = Eq(a*b*c+a*b*d+a*c*d+b*c*d, 8)
h0 = Eq(a*b*c*d, 6720)
hab = Eq(a*b, 80)
# 5 本の連立方程式を a, b, c, d, k について解く
solve([h3, h2, h1, h0, hab], [a, b, c, d, k])

以下が解
4-element Vector{NTuple{5, Sym}}:
 (-10, -8, 7, 12, -178)
 (-10, -8, 12, 7, -178)
 (-8, -10, 7, 12, -178)
 (-8, -10, 12, 7, -178)

いずれの場合も，k = -178 である

Julia での，行列，データフレームの列方向と列方向の統計処理（アクセス）の速度比較

行列，データフレームの列方向と列方向の統計処理（アクセス）の速度比較を行った。
データセットの大きさは 50000行×1000列 である。

列方向の統計処理（アクセス）では，行列もデータフレームも処理時間はほとんど同じである。
行方向の統計処理（アクセス）では，行列では列方向の処理時間の倍になる。データフレームでは，行列に比べて処理時間は400倍になる。

　　　　　行列　　　　　データフレーム
列方向　　0.024035　　 0.024320　　　（秒）
行方向　　0.041314　　16.908246

結論

行方向の統計処理が必要な場合は，データフレームを行列に変換してから分析した方がよい。

using DataFrames, Statistics
using Random; Random.seed!(888)

nr, nc = 50000, 1000;
x = randn(nr, nc);
df = DataFrame(x);

# 列方向
@time colmeansmat = vec(mean(x, dims=1)); # 0.024035 seconds (9 allocations: 8.609 KiB)
@time colmeansdf  = mean.(eachcol(df));   # 0.024320 seconds (1.50 k allocations: 31.672 KiB)
all(colmeansmat == colmeansdf)            # true

# 行方向
@time rowmeansmat = vec(mean(x, dims=2)); # 0.041314 seconds (14 allocations: 391.453 KiB)
@time rowmeansdf  = mean.(eachrow(df));   # 16.908246 seconds (348.64 M allocations: 5.941 GiB, 6.50% gc time)
all(rowmeansmat == rowmeansdf)            # true

行列あるいはデータフレームの行方向の計算

前報は，行列あるいはデータフレームの列方向の計算であった。
今回は，行方向の計算について検討する。

データ規模はやや大きい 50000×1000 を対象とした。

library(microbenchmark)

set.seed(777)
nr = 50000
nc = 1000
x = matrix(runif(nr * nc), nr, nc)
df = as.data.frame(x)

cpu = microbenchmark(
  Reduce.df = Reduce("+", df),
  rowSums.df = rowSums(df),
  apply.df = apply(df, 1, sum),
  for.loop.df = {
      m1 = numeric(nr)
      for (i in 1:nc) {
          m1 = m1 + df[, i]
      }
      m1
  },
  rowSums.mat = rowSums(x),
  apply.mat = apply(x, 1, sum),
  for.loop.mat = {
      M1 = numeric(nr)
      for (i in 1:nc) {
          M1 = M1 + x[, i]
      }
      M1
  },
  unit="ms"
)

par(mar=c(2.7, 6.0, 0.5, 0.5), mgp=c(1.0, 0.5, 0), tcl=-0.2)
plot(cpu,
  col=rep(c("blue", "red"), c(4, 3)),
  las=1, horizontal=TRUE, log = "x",
  xlab="", ylab=""
)
mtext("cpu time in nanoseconds", 1, 1.5)

Unit: milliseconds
         expr        min         lq       mean     median        uq       max neval    cld
    Reduce.df  186.39777  202.03146  243.48758  224.28278  280.4589  455.4433   100  b    
   rowSums.df  335.91715  359.85339  412.49965  409.48675  433.2498  686.4874   100   c   
     apply.df 1804.48749 1971.52982 2147.20970 2092.33050 2304.0389 3154.3865   100      f
  for.loop.df  204.47810  220.42069  274.09268  243.64281  308.5390  510.1865   100  b    
  rowSums.mat   88.33803   89.95369   99.45286   98.05859  108.2877  118.9649   100 a     
    apply.mat 1140.83549 1247.81805 1347.21777 1329.95254 1420.0599 1938.3013   100     e 
 for.loop.mat  387.88910  420.76785  498.32962  498.12949  550.1991  880.5908   100    d  

所見

1. 行列を対象とした rowSums がもっとも速い。for loop は遅いが， apply は相当遅い。
2. データフレームの場合は Reduce と for loop が同程度速く，rowSums は少し遅い。apply はとてつもなく遅い。

なお，行列を対象とする for loop の場合であるが，以下の 2 通りの書き方では M1 = M1 + x[, i] のほうが3倍ほど速いのに注目すべきである。（データフレームの場合も同じ）

理由は，計算中のメモリーアクセスが M1 = M1 + x[, i] は列方向でメモリー上で連続しているのに対し，M1[i] = sum(x[i, ]) は行方向のアクセスで，飛び飛びのメモリーを参照するためである。

system.time({ # 0.277sec.
   M1 = numeric(nr)
      for (i in 1:nc) {
        M1 = M1 + x[, i] # ベクトル演算
      }
})

system.time({ # 0.825sec.
   M1 = numeric(nr)
      for (i in 1:nr) {
        M1[i] = sum(x[i, ]) # 関数の結果を要素に代入
      }
})

system.time({ # 0.140sec.
      m1 = numeric(nr)
      for (i in 1:nc) {
          m1 = m1 + df[, i] # ベクトル演算
      }
})

system.time({ # 0.875sec.
   m1 = numeric(nr)
      for (i in 1:nr) {
        m1[i] = sum(x[i, ]) # 関数の結果を要素に代入
      }
})

9年経ったが，R の「for ループが遅いなんて，誰が言った？」はどうなった？

for ループが遅いなんて，誰が言った？
https://blog.goo.ne.jp/r-de-r/e/49dd48ded7e3605578bbcde95731e412
https://blog.goo.ne.jp/r-de-r/e/4515f9867b017eed1a97cdb8173d8cdc
https://blog.goo.ne.jp/r-de-r/e/1a760a8bee2e3dfb57ffd9235294d4bd
https://blog.goo.ne.jp/r-de-r/e/e1ced6d1474bb9da6a1cf089c15da319
は，2012年の2月末の記事である。
10000×1000のデータフレームでは，forループ≒sapply > colMeans > apply
それを行列に変換したものでは，colMeans > forループ > apply

その後 2018年の3月末に，同じ趣旨で別の人が記事を書いた。

data.frameの高速演算には列ごとならlapply、行ごとならReduceを使おう
https://qiita.com/Atsushi776/items/176426e2195b18eb65b4
行列の場合もデータフレームの場合も，速度は，行列用の関数 > apply > for
データフレームの場合は，lapply が群を抜いて速い。sapply は lapply の結果を simplify するので，ちょっとだけ余分な時間が掛かる。

後者の場合（列ごとの演算の場合）に for は遅いとしているが，

Rのforは遅いと誰が言った？ (data.frameの高速演算には列ごとならlapply、行ごとならReduceを使おうの補足)
https://qiita.com/Atsushi776/items/c31f2345b9c698354c81
行ごとの場合で，100行100,000列 の場合は reduce を使う場合より for ループの方が速い。
ということで，for も場合により速いということ。

その後3年も経って，R も進化したかどうか，どうなっているかをやってみた。
列ごとの集計についてのみ。

> sessionInfo()
R version 4.0.5 Patched (2021-03-31 r80136)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.7

nr = 5000
nc = 100
set.seed(777)
x = matrix(rnorm(nr*nc), nr, nc)
d = as.data.frame(x)

library(microbenchmark)
cpu = microbenchmark(
    lapply.df = lapply(d, sum),
    sapply.df = sapply(d, sum),
    for.loop.df = {
        m1 = numeric(nc)
        for (i in 1:nc) {
            m1[i] = sum(d[,i])
        }
        m1
    },
    colSums.df = colSums(d),
    apply.df = apply(d, 2, sum),
    colSums.mat = colSums(x),
    for.loop.mat = {
        M1 = numeric(nc)
        for (i in 1:nc) {
            M1[i] = sum(x[,i])
        }
        M1
    },
    apply.mat = apply(x, 2, sum),
    unit="ms"
 )

par(mar=c(2.7, 5.8, 0.5, 0.5), mgp=c(1.0, 0.5, 0), tcl=-0.2)
plot(cpu,
 col=rep(c("blue", "red"), c(5, 3)),
    las=1, horizontal=TRUE, log = "x",
    xlab="", ylab=""
)
mtext("cpu time in milliseconds", 1, 1.5)

実行結果

1. 小規模なデータセット

nr = 5000
nc = 100

Unit: milliseconds
         expr       min        lq       mean     median         uq       max neval   cld
    lapply.df  1.047031  1.081145  1.1332266  1.1010855  1.1366530  1.768589   100 a    
    sapply.df  1.074273  1.108415  1.1708123  1.1437340  1.1752400  1.874276   100 a    
  for.loop.df  4.970888  5.328881  6.8997624  5.4409915  5.5839455 65.615935   100   c  
   colSums.df  3.362236  3.518504  3.7168504  3.5814160  3.6827735  9.154761   100  b   
     apply.df 10.747753 12.295381 12.7935086 12.4636335 12.7963960 16.218739   100     e
  colSums.mat  0.420709  0.438020  0.4669311  0.4453005  0.4666935  0.740547   100 a    
 for.loop.mat  7.191930  7.561473  8.6612821  7.7115335  8.0154290 48.668296   100    d 
    apply.mat  7.518857  9.143223  9.4147229  9.3349020  9.5605895 12.260340   100    d 

所見

データフレームの場合は，lapply, sapply が速く，apply は一番遅い。
colSums は以外と遅い。
for loop を使う場合は，colSums より遅い。

行列の場合は，colSums がダントツで速い。
for loop は遅いが，apply よりは速い。

2. 大規模なデータセット

nr = 50000
nc = 1000

Unit: milliseconds
         expr        min        lq       mean     median        uq       max neval   cld
    lapply.df  109.28227  154.7461  238.27326  226.45196  294.1781  606.8323   100 a    
    sapply.df  109.75268  150.9724  223.23908  207.86459  276.3348  641.9150   100 a    
  for.loop.df  126.39996  186.2769  271.16484  245.08860  338.0778 1575.5968   100 a    
   colSums.df  306.50184  440.1110  661.95370  585.96930  825.0703 1725.2760   100  b   
     apply.df 1262.78601 1791.5473 2603.76553 2430.82437 3461.5491 7413.4465   100     e
  colSums.mat   43.77521   64.4854   89.69947   85.70177  108.2785  215.9835   100 a    
 for.loop.mat  464.34522  718.2325 1059.18931  894.83492 1333.1409 3511.8967   100   c  
    apply.mat  960.38814 1426.7567 2049.73518 1778.87375 2604.3161 5885.0181   100    d 

所見

データフレームの場合は，lapply, sapply, for loop が同じくらい速く，apply は一番遅い。
colSums は以外と遅い。

行列の場合は，colSums がダントツで速い。
for loop は遅いが，apply よりは速い。

結論

for はそんなに遅くない。

根拠なき自信

自分ならできる

あの人ならできる

 

でも，なんの根拠もないんだよね

他にもっと強い人がいるということに，目をつぶっているだけ

 

おお，こわいこわい

以て他山の石としよう

SNS を取材ネタにした最近のネットニュースの傾向

タレントなどについて，「○○に対して××という感想がネットで！」などという記事が大変多い。

ほっとけや！どうでもいいじゃん！というのが一般的感想では？

同感しない人はそもそもアンチの投稿はしない。無駄と知っているから。

とくに，何らかの理由でタレントを辞めて，「一般人」になった対象についても，ああだこうだ論評する。一般人なら，ほっとけや！

それらを報じるのが au.one とか livedoor とか。

もう，やめてくれ。そんな記事，もう見たくないんだ。

au.one とか livedoor を表示拒否する訳にはいかんのだ。人質外交だぜ！

分数係数による多項式回帰

分数係数による多項式回帰

格子点を通る多項式を求める場合などで，多項式の係数が分数のままのものが欲しいときがある。
係数が浮動小数点でもよいなら R の lm() を使っても求まるが，特別なプログラムを書いた。

using Statistics, Plots
function RationalLM(x::Vector{Int64}, y::Vector{Int64}, p::Int64 = length(x) - 1)
    n = length(x)
    dat = ones(Rational, n, p + 1);
    for j = 1:p
        dat[:, j] = (x .// big(1)) .^ j
    end
    dat[:, p+1] = big.(y)
    s = cov(dat)
    b = s[1:end-1, 1:end-1] \ s[end, 1:end-1]
    means = mean(dat, dims=1)
    pushfirst!(b, means[end] .- b' * means[1:end-1])
    den = lcm(denominator.(b))
    num = numerator.(b * lcm(den))
    (x=x, y=y, p=p, num=Int64.(num), den=Int64(den), b=b, fb=Float64.(b))
end

predict(obj, x) = [(xi .^ collect(0:obj.p))' * obj.b for xi in x]

function predict(obj)
    pred = predict(obj, obj.x)
    for i = 1:length(obj.x)
        println("x = $(obj.x[i]) \ty = $(obj.y[i]) \tpredict = $(pred[i]) \t= $(Float64(pred[i]))")
    end
end

function summary(obj)
    for i = 1:obj.p + 1
        println("coef$(i-1) \t= $(obj.b[i]) \t= $(obj.num[i]) // $(obj.den) \t= $(obj.fb[i])")
    end
end

function plot_results(obj; width=400, height=300)
    pyplot(grid=false, size=(width, height), label="")
    scatter(obj.x, obj.y, tick_direction=:out)
    minx, maxx = extrema(obj.x)
    margin = (maxx - minx) * 0.1
    x2 = range(minx - margin, maxx + margin, length=1000)
    y2 = predict(obj, x2)
    plot!(x2, y2)
end

使用法

x, y は整数ベクトルであること。

x = collect(0:5);
y = [1,3,7,12,8,15];

あてはめ（多項式の係数を求める）

a = RationalLM(x, y);

多項式の係数を，分数形式，分母を共通とする分数形式，浮動小数点形式でまとめて表示する。

summary(a)
#=
coef0    =    1//1    =   120 // 120    = 1.0
coef1    =  643//60   =  1286 // 120    = 10.716666666666667
coef2    = -151//8    = -2265 // 120    = -18.875
coef3    =  323//24   =  1615 // 120    = 13.458333333333334
coef4    =  -29//8    =  -435 // 120    = -3.625
coef5    =   13//40   =    39 // 120    = 0.325
=#

y の予測値を表示する。

predict(a)
#=
x = 0    y = 1     predict =  1//1     = 1.0)
x = 1    y = 3     predict =  3//1     = 3.0)
x = 2    y = 7     predict =  7//1     = 7.0)
x = 3    y = 12    predict = 12//1     = 12.0)
x = 4    y = 8     predict =  8//1     = 8.0)
x = 5    y = 15    predict = 15//1     = 15.0)
=#

図に示す。

plot_results(a)

第3引数で，多項式の次数を指定できる。

b = RationalLM(x, y, 3);
summary(b)
#=
coef0    =   17//42     =  102 // 252    = 0.40476190476190477
coef1    = 1285//252    = 1285 // 252    = 5.099206349206349
coef2    =   -7//6      = -294 // 252    = -1.1666666666666667
coef3    =    5//36     =   35 // 252    = 0.1388888888888889
=#
predict(b)
#=
x = 0    y = 1    predict =  17//42    = 0.40476190476190477
x = 1    y = 3    predict =  94//21    = 4.476190476190476
x = 2    y = 7    predict = 148//21    = 7.0476190476190474
x = 3    y = 12   predict = 188//21    = 8.952380952380953
x = 4    y = 8    predict = 463//42    = 11.023809523809524
x = 5    y = 15   predict = 296//21    = 14.095238095238095
=#

plot_results(b)

小山田圭吾。辞任！万歳とはいかないぞ！！！

やはり，本人の辞任の前に，解任すべきだったね。

今後とも，著作権料などで悠々自適な生活ができるんだろうけど，今後さらにそれらに水増しされないように，今後は彼の著作物を購入しないとか空桶で彼の歌をリクエストしないとかラジオ・テレビ局にリクエストしないなど，我々もいろいろ対処しないといけないと思う。

まあ，あれこれやっても，優雅に暮らしていけるんだろうな。そんな社会構造なんだろう。幻滅だ。

 

東京五輪組織委、「いじめ自慢」で炎上した小山田圭吾の続投を明言　「倫理観をもって仕事」を評価

も，そうとうな恥っさらしだ

あの小山田圭吾の件

> オンライン署名サイト「Ｃｈａｎｇｅ．ｏｒｇ」で「東京オリパラ開閉会式制作メンバーから小山田圭吾氏の除外を求めます」と題したキャンペーンが１６日にスタート。

というのがあるのだけど，change.org で出てくるページに該当案件のリンクがないように思うんだけど。

教えて！！！

 

受け身じゃダメですね。

以下のようです。

https://www.change.org/p/%E6%9D%B1%E4%BA%AC%E4%BA%94%E8%BC%AA-%E3%83%91%E3%83%A9%E3%83%AA%E3%83%B3%E3%83%94%E3%83%83%E3%82%AF%E7%B5%84%E7%B9%94%E5%A7%94%E5%93%A1%E4%BC%9A-%E6%9D%B1%E4%BA%AC%E3%82%AA%E3%83%AA%E3%83%91%E3%83%A9%E9%96%8B%E9%96%89%E4%BC%9A%E5%BC%8F%E5%88%B6%E4%BD%9C%E3%83%A1%E3%83%B3%E3%83%90%E3%83%BC%E3%81%8B%E3%82%89%E5%B0%8F%E5%B1%B1%E7%94%B0%E5%9C%AD%E5%90%BE%E6%B0%8F%E3%81%AE%E9%99%A4%E5%A4%96%E3%82%92%E6%B1%82%E3%82%81%E3%81%BE%E3%81%99?signed=true

この掲示板，不思議だなあ

なんだか知らないけど，発言直後に12,3のアクセスが必ずある。

一般のかたのアクセスとは思えない。

自動的にサンプルしてるんだろうなあ。本当のアクセス数には反映されない。

小山田圭吾！まだ結論でないの！

言論界も，オリンピック組織委員会も，土日お休みなのか...

結論はわかりきっているのに...

明日一番で，あるいはモーニングショーやお昼の番組でこっぴどく取り上げられないと動かないのかな。

残念な日本に成り下がったものだなあ。

小山田圭吾！歴史に残るゾ！！おめでタイ奴じゃ

オリンピック記念映像（作られるんだろうね）は，未来永劫残る。

開会式の作曲をしたのが誰か，その誰かが学生時代にどんな残虐非道な悪行を働いたのかの記録も残る。

子孫に顔向けできないよね。

それとも，後で編集の結果削除され，なかったことにされる？

それはそれで，淋しいよね？

小山田圭吾。呪われた東京オリンピックにとどめを刺したか？

天皇陛下，菅，橋本，小池，丸川，武藤，政府，自民党，公明党に迷惑を掛けないように，さっさと辞任する方が，よいと思うが？忖度忖度。

小山田圭吾 Keigo Oyamada

謝って済むなら警察はいらん！

英語に訳させたら

If I apologize and finish, I don't need the police.

だって。ばっかじゃない？名訳だなあ。