Julia と R のバージョンの組み合わせで動かないことがある。
現在のところ，以下の組み合わせでは問題ない。
Julia: Version 1.7.3
R: Version 4.2.1 x86_64
M1 用の R だと動かない。
sed を使って単語の置換を行い，結果をファイルに書き出そうと，
run(`sed -e 's/aa/bb/g' foo.txt > out.txt`)
としたら，
LoadError: parsing command `sed -e 's/aa/bb/g' foo.txt > out.txt`:
special characters "#{}()[]<>|&*?~;" must be quoted in commands
とエラーを吐いた。一体どうやって quote したらよいか，ぐぐったら，
run(`sed -e 's/aa/bb/g' foo.txt $(>) out.txt`)
とすればよいみたいだったけど，だめだった。
やれ Homebrew だの pyenve だの conda だの Anaconda だの Miniconda だの
挙句の果てに，
> Pythonを直接インストールする場合、ターミナルを利用せずにインストールできる点はメリットといえるでしょう。
> しかし、それ以外のメリットはありません。
> Homebrew + pyenvや、MinicondaのようにPythonのインストールだけでなく、今後のパッケージ管理の容易さなどを考えると、Pythonの直接インストールはあまりおすすめできません。
> しかし、「とにかくPython3系をインストールして使ってみたい」という方には、簡単にインストールするための手段として有効です。
「それ以外のメリット」ってなんなん。
「パッケージ管理」にあくせくする必要なんかないですよね。
「Pythonの公式サイトにアクセスし、Python3系をダウンロードしてください。」に尽きます。なにをわざわざ苦労する？
数値から文字列への変換
https://qiita.com/dzonesasaki/items/cc33ede78ca679f17f6d
Julia の例が記載されていなかったので...
Python の場合
>>> fVal = 3.141592653589793
>>> str(fVal)
'3.141592653589793'
>>> '{:01.8f}'.format(fVal)
'3.14159265'
>>> f'{fVal:01.8f}'
'3.14159265'
Julia では
julia> fVal = pi
π = 3.1415926535897...
julia> string(fVal)
"π"
ulia> string(1pi)
"3.141592653589793"
julia> using Printf
julia> @printf("%.8F", fVal)
3.14159265
https://mac.r-project.org/
もう，2ヶ月近くか，Status OK でなかったのだけど，今日見たら，（4.3.0 も含め）x86_64 用も arm64 用も OK になっていた。
まあ，説明するのも面倒なのですが，学校で習った四捨五入は「不適切な方法」なので，忘れましょう。
JIS Z 8401（数値の丸め方）（PDF、75KB）
https://www.erca.go.jp/fukakin/before/pdf/suuchi.pdf
あなたが思ったのと違った結果になっても，それは，「あなた（学校教育あるいは文部科学省）が間違っている（不適切だ）から」なのです（残念！！！）。
R の round の example
> example(round)
round> round(.5 + -2:4) # IEEE / IEC rounding: -2 0 0 2 2 4 4
[1] -2 0 0 2 2 4 4
round> ## (this is *good* behaviour -- do *NOT* report it as bug !)
round>
round> ( x1 <- seq(-2, 4, by = .5) )
[1] -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
round> round(x1) #-- IEEE / IEC rounding !
[1] -2 -2 -1 0 0 0 1 2 2 2 3 4 4
round> x1[trunc(x1) != floor(x1)]
[1] -1.5 -0.5
round> x1[round(x1) != floor(x1 + .5)]
[1] -1.5 0.5 2.5
round> (non.int <- ceiling(x1) != floor(x1))
[1] FALSE TRUE FALSE TRUE FALSE TRUE FALSE TRUE FALSE TRUE FALSE TRUE FALSE
round> x2 <- pi * 100^(-1:3)
round> round(x2, 3)
[1] 0.031 3.142 314.159 31415.927 3141592.654
round> signif(x2, 3)
[1] 3.14e-02 3.14e+00 3.14e+02 3.14e+04 3.14e+06
Python で同じことを書いてみる。
>>> np.around(0.5 + np.arange(-2, 5))
array([-2., -0., 0., 2., 2., 4., 4.])
>>> x1 = np.arange(-2, 4.5, step=0.5)
>>> x1
array([-2. , -1.5, -1. , -0.5, 0. , 0.5, 1. , 1.5, 2. , 2.5, 3. ,
3.5, 4. ])
>>> np.around(x1)
array([-2., -2., -1., -0., 0., 0., 1., 2., 2., 2., 3., 4., 4.])
>>> x1[np.trunc(x1) != np.floor(x1)]
array([-1.5, -0.5])
>>> x1[np.around(x1) != np.floor(x1 + 0.5)]
array([-1.5, 0.5, 2.5])
>>> non_int = np.ceil(x1) != np.floor(x1)
>>> non_int
array([False, True, False, True, False, True, False, True, False,
True, False, True, False])
>>> x2 = np.pi * 100.0**(np.arange(-1, 4))
>>> np.around(x2, 3)
array([3.10000000e-02, 3.14200000e+00, 3.14159000e+02, 3.14159270e+04,
3.14159265e+06])
>>> print(np.around(x2, 3))
[3.10000000e-02 3.14200000e+00 3.14159000e+02 3.14159270e+04
3.14159265e+06]
>>> print(*np.around(x2, 3))
0.031 3.142 314.159 31415.927 3141592.654
R の signif を Python の関数として書く。
>>> def signif(x, dig):
... return np.around(x, dig - int(np.ceil(np.log10(abs(x)))))
R のようには，いわゆる「科学的浮動小数点表示」にはならない。
>>> result = [signif(x, 3) for x in x2]
>>> print(*result)
0.0314 3.14 314.0 31400.0 3140000.0
Julia で同じことを書いてみる。
結果を整数 Int にしないと -0.0 が出現する。
round.(0.5 .+ -2:4) |> println # [-2.0, -0.0, 0.0, 2.0, 2.0, 4.0]
round.(Int, 0.5 .+ -2:4) |> println # [-2, 0, 0, 2, 2, 4]
x1 = -2:0.5:4
round.(x1) |> println # [-2.0, -2.0, -1.0, -0.0, 0.0, 0.0, 1.0, 2.0, 2.0, 2.0, 3.0, 4.0, 4.0]
x1[trunc.(x1) .!= floor.(x1)] |> println # [-1.5, -0.5]
x1[round.(x1) .!= floor.(x1 .+ 0.5)] |> println # [-1.5, 0.5, 2.5]
non_int = ceil.(x1) .!= floor.(x1)
non_int |> println # Bool[0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0]
x2 = π .* 100.0 .^(-1:3)
round.(x2, digits=3) |> println # [0.031, 3.142, 314.159, 31415.927, 3.141592654e6]
Julia にも R の signif 関数を定義する。
signif(x, dig) = round(x, digits=dig - ceil(Int, log10(abs(x))))
最後の数字のように，巨大数になると「科学的浮動小数点表示」になるが，これは print 関数の働きの結果であろう。
signif.(x2, 3) |> println # [0.0314, 3.14, 314.0, 31400.0, 3.14e6]
Python 公式ドキュメントに
> 注釈 浮動小数点数に対する round() の振る舞いは意外なものかもしれません: 例えば、 round(2.675, 2) は予想通りの 2.68 ではなく 2.67 を与えます。これはバグではありません: これはほとんどの小数が浮動小数点数で正確に表せないことの結果です。詳しくは 浮動小数点演算、その問題と制限 を参照してください。
と書いてあるそうだ。この説明は正確ではない。
2.5 も 3.5 も「浮動小数点数で正確に表せる」が，round の結果は異なる。
>>> round(2.5)
2
>>> round(3.5)
4
2.5 は 16 進表現で (2.8)16 で，10 進数に変換するのは 2 * 16^0 + 8 * 16^(-1) = 2 + 0.5 である。
3.5 は 16 進表現で (3.8)16 で，10 進数に変換するのは 3 * 16^0 + 8 * 16^(-1) = 2 + 0.5 である。
つまり，小数部が 1/2^1, 1/2^2, 1/2^3, 1/2^4, 1/2^5, 1/2^6... = 0.5, 0.25, 0.125, 0.0625, 0.03125, 0.015625... の任意の組み合わせで表せる場合は，「浮動小数点数で正確に表せる」。
0.1640625 = 1/2^3 + 1/2^5 + 1/2^7 であり，16進数で表せば (0.2A)16，2進数で表せば (0.0010101)2 つまり，2進数の小数点以下 3, 5, 7 が 1 の数であり，正確に表すことができるのである。
Python の公式文書は（翻訳者がミスしただけかもしれないが）ダメダメであるが，他の言語ではどうだろうか?
R: `? round` において
‘Details’ about “round to even” when rounding off a 5
Rounding to decimal digits in binary arithmetic is non-trivial (when
digits != 0) and may be surprising. Be aware that most decimal fractions are not exactly representable in binary double precision. In R 4.0.0, the algorithm for round(x, d), for d>0d>0, has been improved to measure and round “to nearest even”, contrary to earlier versions of R (or also to sprintf() or format() based rounding).
と書いているように，「四捨五入をしようとしたとき，結果（の末尾）が偶数になるように丸める」のである。
つまり，2.5 を「学校で習った四捨五入」をすると 3 になるのであるが，これは奇数なので 2 に丸める（2.5 は 2 と 3 のちょうど真ん中なので，2 に丸めても 3 に丸めてもよいのだけど結果が偶数になるように 2 に丸めるということ）。
3.5 を「学校で習った四捨五入」をすると 4 になるのであるが，これは偶数なので 4 のままにする。
Julia でも，being rounded to the nearest even integer と書いてある。
The RoundingMode r controls the direction of the rounding; the default is RoundNearest, which rounds to the nearest integer, with ties (fractional values of 0.5) being rounded to the nearest even integer. Note that round may give incorrect results if the global rounding mode is changed (see rounding).
こうなると，Python はどうなんだ?と思う。help(np.around) で表示される文書に，ちゃんと書いてある（np.round より np.around をつかえといわれる）。
For values exactly halfway between rounded decimal values, NumPy rounds to the nearest even value. Thus 1.5 and 2.5 round to 2.0, -0.5 and 0.5 round to 0.0, etc.
これに続いて更に詳しく書かれている。ぜひちゃんと読むべし。
冒頭に掲げた，寝言（の原文）はどこに書いてあるのかな?
「正規序数値の期待値」とは聞き慣れない訳語かもしれないが，Expected values of normal order statistics のこと。
H. Leon Harter（1961）Expected Values of Normal Order Statistics, Biometrika, 48, 151-165.
https://www.jstor.org/stable/2333139
大本に近いのは以下
Algorithm AS 177: Expected Normal Order Statistics (Exact and Approximate)
http://imaging.mrc-cbu.cam.ac.uk/statswiki/FAQ/Corder?action=AttachFile&do=get&target=royston.pdf
FORTRAN 90 に書き換えたのが
as177.f90
https://jblevins.org/mirror/amiller/as177.f90
であるが，correc 関数中の定数 C3(4) の 2.157E6 が 2.1576E6 と誤転写されている模様（C に書き換えられたものもあったが 2.157E6 だった）。
なにはともあれ，Julia で最適化して（?）以下のコードを得る。
nscor1 が正確なもの。nscore2 は近似解。
using SpecialFunctions
using Distributions
function nscor1(n)
"""
Exact calculation of Normal Scores
"""
NSTEP = 721
work = zeros(4, NSTEP)
init!(work)
n2 = div(n, 2)
s = zeros(n2)
h = 0.025
c1 = logfactorial(n)
d = c1 - log(n)
for i in 1:n2
i1 = i-1
ni = n-i
c = c1 - d
scor = 0
for j in 1:NSTEP
scor += exp(work[2, j] + i1 * work[3, j] + ni * work[4, j] + c) * work[1, j]
end
s[i] = scor * h
d += log(i / ni)
end
return s
end
function init!(work)
xstart = -9
h = 0.025
pi2 = -0.918938533
xx = xstart
for i in 1:size(work, 2)
work[1, i] = xx
work[2, i] = pi2 - 0.5 * xx^2
work[3, i] = log(ccdf(Normal(), xx))
work[4, i] = log(cdf(Normal(), xx))
xx = xstart + i * h
end
end
近似解 nscor2 は以下
function nscor2(n)
"""
Calculates approximate expected values of normal order statistics.
claimed accuracy is 0.0001, though usually accurate to 5-6 dec.
"""
eps = [0.419885, 0.450536, 0.456936, 0.468488]
dl1 = [0.112063, 0.121770, 0.239299, 0.215159]
dl2 = [0.080122, 0.111348, -0.211867, -0.115049]
gam = [0.474798, 0.469051, 0.208597, 0.259784]
lam = [0.282765, 0.304856, 0.407708, 0.414093]
bb = -0.283833
d = -0.106136
b1 = 0.5641896
n2 = div(n, 2)
s = zeros(n2)
s[1] = b1
n == 2 && return
k = 3
n2 < k && (k = n2)
for i in 1:k
e1 = (i - eps[i]) / (n + gam[i])
e2 = e1^lam[i]
s[i] = e1 + e2 * (dl1[i] + e2 * dl2[i]) / n - correct(i, n)
end
if n2 != k
for i in 4:n2
l1 = lam[4] * bb / (i + d)
e1 = (i - eps[4]) / (n + gam[4])
e2 = e1^l1
s[i] = e1 + e2 * (dl1[4] + e2 * dl2[4]) / n - correct(i, n)
end
end
for i in 1:n2
s[i] = -quantile(Normal(), s[i]) # -ppnd(s[i])
end
return s
end
function correct(i, n)
"""
Calculates correction for tail area of the i-th largest of n order statistics.
"""
c1 = [9.5, 28.7, 1.9, 0.0, -7.0, -6.2, -1.6]
c2 = [-6195.0, -9569.0, -6728.0, -17614.0, -8278.0, -3570.0, 1075.0]
c3 = [9.338e4, 1.7516e5, 4.1040e5, 2.157e6, 2.376e6, 2.065e6, 2.065e6]
mic = 1e-6
c14 = 1.9e-5
i * n == 4 && return c14
(i < 1 || i > 7) && return 0
(i != 4 && n > 20) && return 0
(i == 4 && n > 40) && return 0
an = 1 / n^2
return (c1[i] + an * (c2[i] + an * c3[i])) * mic
end
呼び出し方はいずれも同じ n だけを引数にして呼ぶ。
nscore1(30)
15-element Vector{Float64}:
2.0427608445896133
1.6155999007846518
1.3648120475550987
1.1785538961988162
1.0260849654763988
0.8943875579926721
0.7766582371097811
0.6688521087544627
0.5683389891912317
0.47328800401844007
0.38235131764319513
0.29448658895879626
0.2088497011984587
0.12472539951340163
0.04147914900072655
nscor2(30)
15-element Vector{Float64}:
2.042747144755947
1.6155640427483784
1.3648321622142314
1.1749378349295097
1.023401473812151
0.892300080600137
0.7750026028133499
0.667520049660839
0.5672570676344856
0.47240515390672405
0.38163143261223
0.29390400347499945
0.20838637645235597
0.12436915904576364
0.04122259874887312
「2020年1月から2022年5月まで、連続の月次データを生成する」ということで，持って回ったプログラムが紹介されている。
library(tidyverse)
library(lubridate)
year_month <- seq(as.Date("2020-01-01"), as.Date("2022-05-01"),by="1 day") %>% floor_date(unit="months") %>% unique()
seq(as.Date("2020-01-01"), as.Date("2022-05-01"), by="1 month")
で十分だし，
seq(as.Date("2020-01-01"), as.Date("2022-05-01"), by="1 week")
seq(as.Date("2020-01-01"), as.Date("2022-05-01"), by="1 quarter")
seq(as.Date("2020-01-01"), as.Date("2022-05-01"), by="2 week)
seq(as.Date("2020-01-01"), as.Date("2022-05-01"), by="5 day")
などなど，なんでもありなんだけどなあ
today <- Sys.Date() seq(today, length.out=10, by="1 week")
なんてこともできる。
「オンラインヘルプを見よ！」と言われても，? seq じゃあ出て来ない。? Date でヒントが得られる。まあ，そういうもんだろうけど。
PVアクセスランキング にほんブログ村
