R の fivenum 関数を Julia に移植する。 function fivenum!(x) sort!(x) n = length(x) n4 = floor(Int, (n + 3) / 2) / 2 d = [1, n4, (n + 1) / 2, n + 1 - n4, n] @. 0.5 * (x[floor(Int, d)] + x[ceil(Int, d)])end fivenum! は，最小値，下ヒンジ，中…

Julia では π は特殊な型を持つ typeof(π) # Irrational{:π} なので， π # π = 3.1415926535897... という表示になる。 ... に意味がある。 big(π) あるいは BigFloat(π) は 3.141592653589793238462643383279502884197169399375105820974944592307816406286…

どうという計算ではないのに，結果が違う。 実数演算にはご注意ということではあるが。 x = 0.0:0.1:1.0 n が 123 のときは，両辺等しい。 n = 123n*x .== n*collect(x) 11-element BitVector: 1 1 1 1 1 1 1 1 1 1 1 n が 122 のときは，等しくない n = 122…

また，最近，「妖怪ダイナマイトプロット」 dynamite plot, plunger plot, dynamite plumger が出現か これを勧める，描画するプログラムを掲出するということがあるので，再度，警告を発す。 A/B テストの周辺（中心？）で，御託を並べるダメ出し：ダイナマ…

いわゆる電卓サイトで 4π÷2π を計算するとどうなるか。実際にやってみた。 Julia だと 4π / 2π ==> 2 なのだが， 4π / 2π == 4 * π / 2 * π ==> 19.739208802178716 とするものや，そもそも 4π という入力を許さないものなど，さまざまだ。 まあ，多くのプロ…

missing は他の言語の NaN ではない。NaN も他の言語の NaN ではない？ b = NaNisnan(b) || ismissing(b) # trueisnan(b) | ismissing(b) # true c = missingisnan(c) || ismissing(c) # TypeError: non-boolean (Missing) used in boolean contextisnan(c) …

コイントスの結果は2項分布になるので，B(n, p-0.5) の n が大きくなれば正規分布に近づくというのは前に述べた。 コインの裏表を 1/0 で表して合計すると（平均値を求めても同じであるが）考えれば，一様分布に従う n 個のデータの和ということで，中心極限…

自由度が∞の場合は，正規分布である。 自由度が1の場合はコーシー分布で，見た目からは信じられないだろうが，平均値は定義されない。 using Plots, PlotThemes, Rmath function tdistribution(;fps=10) x = -4:0.1:4 theme(:gruvbox_light) pyplot(grid=fal…

自由度が 1, 2 のときは形がちょっと違う。 作成するプログラムは以下の通り。 using Plots, PlotThemes, Rmath function chisquaredistribution(;fps=2) x = 0:0.1:100 theme(:orange) pyplot(grid=false, label="", size=(400, 300)) anim = @animate for …

収穫物：前の記事で，2つのグラフをまとめて gif ファイルにする方法がわからないとか，愚痴っていたけど，簡単だった。自己解決。 plt1 = foo1();plt1 = bar1();plt2 = baz2()plt2 = boo()plot(plt1, plt2, layout=(n, m)) とすればよいだけだった。 という…

ポアソン分布においては，ポアソン定数 λが大きくなるにつれ，正規分布に近づく。 ただし，二項分布が正規分布に近づくのに比べると，遅い。 最終時点での結果は以下の通り。正規分布曲線にはまだ十分近づいているとは言いがたい。 using Plots, PlotThemes,…

二項分布 B(n, p) において，n がそんなに大きくない場合は，p が 0.5 でないときには分布はかなり歪んでいる。しかし，n が大きくなると，だんだん正規分布に近づいていく。 その様子を示したのが以下の GIF アニメーション。最後の状態 n = 75 のときの場…

散布図を見ただけで，ピアソンの積率相関係数のだいたいの値をいえるようになると，経験値が上がるかもしれない... 作成する Julia プログラムは以下の通り。 using LinearAlgebra, StatsBase, Statisticsusing Plots, StatsPlots, PlotThemesusing Random f…

R には動画画像ファイルの作成のために animation ライブラリがあるが，Julia にも動画 gif ファイルを作成する機能がある。 試しに，モンテカルロ法で π を求める動画を作ってみた。 以下は，2つの gif ファイルを左右に並べて表示しているので（環境によっ…

自由落下における速度，距離，時間 重力加速度(g) ≒ 9.8 m/s^2 時間(t;s)と速度(v;m/s) 次元解析 (m/s) = (m/s^2) * （s) v = gt 時間(t;s)と落下距離(h;m) 次元解析 (m) = (m/s^2) * （s^2) h = gt^2 / 2 落下距離(h;m)に必要な時間(t;s) 次元解析 （s) = √…

Mac OS Catalina, バージョン 10.15.7 で gmp 0.6−2 が動かない。爆弾（エラー）が出る。 いろいろ探して， M1チップ搭載のMacだとRのgmpパッケージが動かない？https://rion778.hatenablog.com/entry/2021/03/02/214104 に， install.packages("gmp", type …

図形問題の解 鶏を割くに牛刀を以てす ではあるが， 件の図で，f を原点(0, 0) として，fb を f = 5/2 *xfd を g = x / 4bd を h = -2x + 9とし，f - g, h - g の定積分を求め，2倍する Julia の SymPy でやってみよう using SymPy @syms f g h x f = 5/2 * …

https://www.youtube.com/watch?v=51Ru7F2uMms&feature=emb_rel_end 小学生でも解けるかな？ っていうか，私が小学生並みなだけ。 下図のように頂点の記号を定める。平行四辺形の2つの頂点 b, d を結ぶ。ab の延長線とed の延長線の交点が c。af, ce は垂直…

R の bincombinations() R の bincombinations(p) は，0/1 からなる 2^p 個の長さ p のベクトルを生成する。 Julia に移植すると以下のようになる。 function bincombinations(p) ary = Array{Int,2}(undef, 2^p, p); for n = 1:p ary[:, n] = repeat(vcat(f…

(x - a)(x - b)…(x - z) を展開するとどうなるかということで，答えは 0，なぜならば途中に (x-x) すなわち 0 があるでしょうということで，まあ，トンチ問題，引っかけ問題ではあるのだけど。 しかし，これを Julia の SymPy で解こうとすると，「え！！」…

前報の無理数の連分数近似というのは，それ自体はあまりありがたいものでもない。 しかし，ペル方程式 Pell's equation を解くときに力を発揮する。 ペル方程式とは， n を平方数ではない自然数として、未知整数 x, y について x^2 - n*y^2 = 1 の形のディオ…

連分数展開とは，たとえば x の平方根を以下のような連続する分数の形式で表現することである。 √x = a0 + 1/(a1 + 1/(a2 + 1/(a3 + 1/(a4 + 1 / (a5 + ε)))))εはさらに 1/(a6 + 1/(a7 + 1/...)) と同じように無限に繰り返されるが，だんだんと近似に寄与す…

qiita.com/tags/python で，i += 1 より i = i+1 のほうが速いというのがあったが，そもそも，10000000 回！！もやって所要時間が 0.69685sec と 0.67206sec なんだから，誤差範囲でしょう。 ちなみに，Julia でやったら，誤差範囲で同じようなものだった。…

Julia の Primes モジュールは，R や Python で gmp パッケージでサポートされているものと同じだ。 using Primes 素因数分解 factor(n::Integer), factor(ContainerType, n::Integer) 結果の表示法 特に指定しない場合は，数式として見やすい結果が得られる…

同じ文字列が何個続くか勘定する。 aabbbcccc なら，a が 2個，b が 3個，c が 4個なので，([2, 3, 4], ["a", "b", "c"]) のようなタプルを返す。 function rle(x) """ Run Length Encoding Compute the lengths and values of runs of equal values in a ve…

using Polynomialsusing Plotsgr(xlabel="\$x\$", ylabel="\$y\$", size=(400, 300)) 多項式の定義，多項式の値，多項式の解 Polynomial() の目的は，1つには与えた係数ベクトルから「多項式 Polynomial{Float64, :x}」を作ることである。係数ベクトルは定数…

x = range(0, 5, length=6)y = @. exp(-x)#=6-element Vector{Float64}: 1.0 0.36787944117144233 0.1353352832366127 0.049787068367863944 0.01831563888873418 0.006737946999085467 =# というデータがある。y を x の多項式 a0 + a1*x + a2*x^2 + ... で…

#==========Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。以下のプログラムを Julia に翻訳してみる。 陽性反応適中率の差の検定http://aoki2.si.gunma-u.ac.jp/JavaScript/ppv.html ファイル名: ppv.jl 関数名: ppv 翻訳…

いま，たまたま 「ggplot2 って，ちっとも簡単じゃないよ！」が，参照されたようだが，ggplot2 で描画される画像と，VegaLite の画像を比べてみた。 VegaLite は巨大なパッケージのようで，指定もよく分からない。 試行錯誤の結果， グリッドを消すのは，axi…