裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

Julia の iterator

2021年01月18日 | ブログラミング

Python でも使われる enumerate(), zip() の他にもいろいろある。

julia> A = 1:5
julia> B = ["one", "two", "three", "four", "five"]

タプルで取り出す場合は,変数 1 個を使う。取り出した後はタプルとして使ってもよいし,ここの変数へバラしてから使ってもよい。

julia> for c in zip(A, B)
           a, b = c
           println("$c, $a, $b")
       end
(1, "one"), 1, one

(2, "two"), 2, two
(3, "three"), 3, three
(4, "four"), 4, four
(5, "five"), 5, five

各要素へ取り込む場合は,Python と異なり,(  ) でくくらないとエラーになる。

julia> for (a, b) in zip(A, B)
           println("$a, $b")
       end
1, one
2, two
3, three
4, four
5, five

enumerate() も,取り出し変数は (  ) でくくること。

julia> for (index, value) in enumerate(B)
           println("$index, $value")
       end
1, one
2, two
3, three
4, four
5, five

面白そうな iterator を一つあげておく。汎用の outer 関数を作るときに使えるかな??

julia> a = collect(Iterators.product(1:2, 3:5))
2×3 Array{Tuple{Int64,Int64},2}:
 (1, 3)  (1, 4)  (1, 5)
 (2, 3)  (2, 4)  (2, 5)

julia> println(a)
[(1, 3) (1, 4) (1, 5); (2, 3) (2, 4) (2, 5)]

julia> for x in a
           println(x)
       end
(1, 3)

(2, 3)
(1, 4)
(2, 4)
(1, 5)
(2, 5)

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia で微分,求解

2021年01月18日 | ブログラミング

某所で,sin(θ) + sqrt(3) * cos(θ) の最大値を求めよ。但し 0 ≦θ≦π/2。を解いていたのだが,中途半端だったので。

julia> using SymPy

julia> @syms θ f
(θ, f)

julia> f = sin(θ) + sqrt(3) * cos(θ)
sin(θ) + 1.73205080756888⋅cos(θ)

julia> x = solve(diff(f)) # 微分して 0 になるときの x の値を求める
1-element Array{Sym,1}:
 0.523598775598298 # π/6

julia> f.subs(θ, x[1]) # 元の式に代入すればその時の値(最大値)が求まる。
2.00000000000000

julia> simplify(diff(f))
-1.73205080756888⋅sin(θ) + cos(θ)

julia> using Plots # 図を描いて確認する

julia> g(x) = sin(x)+sqrt(3)*cos(x) # 元の式
g (generic function with 1 method)

julia> h(x) = -sqrt(3)*sin(x)+cos(x) # 微分した式
h (generic function with 1 method)

julia> xax = 0:0.01:pi/2
0.0:0.01:1.57

julia> plot(xax, g.(xax), label="sin(θ) + sqrt(3) * cos(θ)") # 元の式(水色)

julia> plot!(xax, h.(xax), label="-sqrt(3)*sin(θ) + cos(θ)") # 微分した式(赤)

julia> vline!([x], lw=0.3, label="") # 微分した式が 0, 元の式が最大値になるときの x を通る垂直線

julia> hline!([0.0], lw=0.3, label="") # x 軸

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia で微分

2021年01月18日 | ブログラミング

うう〜ん。これって,Python の sympy 使ってるだけだよな〜〜〜

θによる媒介変数表示

x = sin θ
y = -log tan(θ/2) - cos θ   (0 < θ < π/2)

θ=π/3 における d2y / dx2 の値を求めよ。log は自然対数関数。

using SymPy

@syms a t u x y z
x = sin(t)
y = -log(tan(t/2)) - cos(t)
z = diff(diff(y)/diff(x))/diff(x)
u = simplify(z)
a = subs(u, t, pi/3)
round(Int, N(10a - a)) // 9

答え 8//3

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia でスピアマンの順位相関係数行列

2021年01月17日 | ブログラミング

ご存じのこととは思いますが,念の為

Julia で 2 変数間のスピアマンの順位相関係数を求めるには,

using StatsBase, Statistics
corspearman([1,2,3,4,5], [3,2,1,2,3])

です。

スピアマンの順位相関係数行列を求めるとき,

function spearmancor0(x::Array{Float64,2})
    nc = size(x, 2)
    r = ones(nc, nc)
    for i = 2:nc
        for j = 1:i
            r[i, j] = r[j, i] = corspearman(x[:, i], x[:, j])
        end
    end
    return r
end

としてはいけません

corspearman は データ行列も受け付けます。
以下のように定義した関数がやるように,データ行列を列ごとに平均順位を付け,そのあとで cor 関数を呼びます。

function spearmancor(x::Array{Float64,2})
    rank = similar(x, Float64)
    for i = 1:size(x, 2)
        rank[:, i] = tiedrank(x[:, i])
    end
    return cor(rank)
end

のようにしましょう。

時間計測すると,以下のようになります。

x = randn(10000, 100);
@time b = spearmancor0(x); # 8.367428 seconds (70.69 k allocations: 2.260 GiB, ...
@time a = spearmancor(x);  # 0.097368 seconds (714 allocations: 38.250 MiB, ...
@time c = corspearman(x); # 0.094456 seconds (1.06 k allocations: 38.337 MiB
a == b # true
a == c # true

90 倍ほど速いです。当然ですが,得られる相関係数行列は同じです。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

セイウチ演算子(:=)って今更...(その2)

2021年01月17日 | ブログラミング

いやあ,うっかりしてたわ。

R でセイウチ演算子での代入箇所に複文が書けるって書いていたのに,Julia のこと書かなかったわ。

Julia では,複文は begin 〜 end で囲むので次の様にして,複文を書くことができる。

どんなに複雑な式を何個書いても問題ない。

が,もっとも大事な注意があることを書いてなかったわ。

a = [1,2,3,4,5]
if begin n = length(a); w=20*pi end > 10
    println("List is too long ($n elements, expected <= 10)")
else
    println("variable n is $n, w = $w")
end
println(n, ", ", w)

List is too long (5 elements, expected <= 10), w = 62.83185307179586
5, 62.83185307179586

おや,5 elements, expected <= 10 になってますよ。

あわてちゃいけない。複文の評価は,最後の文の評価結果なの!!つまり, w > 10 をやったことになってるの。

位置を考えるか,最後に n を置くか。以下のようにする。

if begin n = length(a); w=20*pi; n end > 10

==========

R も同じ。n > 10 をやりたいなら,複文の最後に n を置く

a = c(1,2,3,4,5)
if ({n = length(a); w = 20*pi; n} > 10) {
    cat(sprintf("List is too long (%d elements, expected <= 10), w = %f\n", n, w))
} else {
    cat(sprintf("variable n is %d, w = %f\n", n, w))
}
cat(n, w)

variable n is 5, w = 62.831853
> cat(n, w)
5 62.83185

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

セイウチ演算子(:=)って今更...

2021年01月17日 | ブログラミング

Python 3.8 になって,使えるようになって感涙にむせんでいる人もいるようなんだけど...

a = [1,2,3,4,5]
if (n := len(a)) > 10:
    print(f"List is too long ({n} elements, expected <= 10)")
else:
    print(f"variable n is {n}")
print(n)

実行結果

variable n is 5
5

しかし,これと同じ機能は多くの言語では古くから使えていた
というか,使えていたので,使い方を間違って問題を起こすこともあった。

R では,:= を使うときにも = を使えば良い。ただし,'(代入文)' のように '( )' でくくることを忘れてはいけない。

a = c(1,2,3,4,5)
if ((n = length(a)) > 10) {
    cat(sprintf("List is too long (%d elements, expected <= 10)\n", n))
} else {
    cat(sprintf("variable n is %d\n", n))
}
print(n)

実行結果

variable n is 5
> print(n)
[1] 5

R では if 文の条件式を ( ) でくくらなければならないので,Python と同じようには見えないかも知れないが,外側の( ) を取り除けば Python のセイウチ演算子と同じなのだ。

Julia だともっとわかりやすいだろう(条件式を ( ) でくくる必要がない)。
Python のセイウチ演算子の位置には,代入演算子 = を使っているだけだ。
セイウチ演算子なんてなくても,ちっとも困らない。

a = [1,2,3,4,5]
if (n = length(a)) > 10
    println("List is too long ($n elements, expected <= 10)")
else
    println("variable n is $n")
end
println(n)

実行結果

variable n is 5
5

古いというほど古くはない C でも以下のようになる。
代入式を(n = sizeof(a)/sizeof(int)) のように ( ) でくくるだけ。

#include
#include
int main(void) {
      int n;
      int a[5] = {1,2,3,4,5};
      if ((n = sizeof(a)/sizeof(int)) > 10) {
            printf("List is too long (%d elements, expected <= 10)\n", n);
      }
      else {
            printf("variable n is %d\n", n);
      }
      printf("n = %d\n", n);
      return EXIT_SUCCESS;
}

実行結果

variable n is 5
n = 5

ところで,R はもっと進んでいるんだ。

a = c(1,2,3,4,5)
if ({n = length(a)} > 10) {
    cat(sprintf("List is too long (%d elements, expected <= 10)\n", n))
} else {
    cat(sprintf("variable n is %d\n", n))
}
print(n)

のように,代入文を { } でくくっても良い
それは,この位置に複文を書けるということだ。

以下のようなこともできる。すきなだけ複数の複雑な式を書ける。

a = c(1,2,3,4,5)
if ({n = length(a); m = 2*pi} > 10) {
    cat(sprintf("List is too long (%d elements, expected <= 10)\n", n))
} else {
    cat(sprintf("variable n is %d\n", n))
    cat(sprintf("2*pi = %f\n", m))
}
cat(n, m)

実行結果

variable n is 5
2*pi = 6.283185
> cat(n, m)
5 6.283185

しかし,同じように { } で複文が書ける C では,この位置に複文を書くとエラーになる。

もっと古い AWK でも同じように1つの文は書けるが,やはり複文は書けない。

BEGIN {
      a = 2
      b = 6
      if ((n = a * b) > 0) {
            print "true", n
      }
      else {
            print "false", n
      }
      print "outer", n
}

実行結果

true 12
outer 12

要するに,新しい演算子を作るかどうかということではなく,文法がどうか(記述されたものをどのように解釈するか)ということになるだろう。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

数学の問題を Julia で解く

2021年01月17日 | ブログラミング

某所に Python による解答例があったので,Julia でやるならどうなるか書いてみた。

Q. 1. 2018n ≡ 2 mod(1000) を満たす正の整数nの最小値を求めよ。

for n = 1:99999
  if mod(2018n, 1000) == 2
    println(n) # 389
    break
  end
end

Q. 2. tan(2 Arctan (1/3)+Arctan(1/12))を求めよ。

println(tan(2atan(1//3)+atan(1//12))) # 0.8888888888888888

a = tan(2atan(1//3)+atan(1//12))
println(Int(10a - a) // 9) # 8//9

Q. 3. xyz空間内の4点(1,-4,1),(2,2,2),(2,-6,-3),(3,-2,-1)を頂点に持つ四面体の体積を求めよ。

sarasu(x1, y1, z1, x2, y2, z2, x3, y3, z3) =
    abs((y1*z2*x3+z1*x2*y3+x1*y2*z3-z1*y2*x3-x1*z2*y3-y1*x2*z3) / 6)

a1 = [1, -4, 1];
b1 = [2, 2, 2];
c1 = [2, -6, -3];
d1 = [3, -2, -1];
x1, y1, z1, x2, y2, z2, x3, y3, z3 = vcat(b1 .- a1, c1 .- a1, d1 .- a1);
println(abs((y1*z2*x3+z1*x2*y3+x1*y2*z3-z1*y2*x3-x1*z2*y3-y1*x2*z3) / 6)) # 3.0
println(sarasu(x1, y1, z1, x2, y2, z2, x3, y3, z3)) # 3.0

別解

A = hcat(b1 .- a1, c1 .- a1, d1 .- a1);
using LinearAlgebra
println(abs(det(A)) / 6) # 3.0

Q. 4.1 3次正方行列 A=[4 1 1; 1 2 1; 1 1 2] について固有値を求めよ。

A = [4 1 1; 1 2 1; 1 1 2];
using LinearAlgebra
println(round.(Int, eigvals(A))) # [1, 2, 5]

Q. 4.2 E=[1 0 0; 0 1 0; 0 0 1] のとき,A^3-8A^2+18A-12E を求めよ

E = Array{Int}(I, 3, 3);
println(A^3 .- 8A^2 .+ 18A .- 12E) # [2 1 1; 1 0 1; 1 1 0]

Q. 5 2つの機械A, Bが、1年間に故障する回数は,平均がそれぞれ2.5回、4.5回のポアソン分布に従う。A,Bの故障は独立である。
Q. 5.1 機械Aが1年間に4回以上故障する確率を求めよ。

using Distributions
println(round(ccdf(Poisson(2.5), 3), digits=3)) # 0.242

Q. 5.2 機械A,Bが合わせて8回以上故障する確率を求めよ。

A = Poisson(2.5)
B = Poisson(4.5)
answer = 1.0
for i = 0:7
    for j = 0:i
        answer -= pdf(A, i-j) * pdf(B, j)
    end
end
println(round(answer, digits=3)) # 0.401

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia で一元配置分散分析(データフレーム)

2021年01月15日 | ブログラミング

Julia でANOVA」 という記事がある。記事の日付は2019年3月15日とわずか2年足らず前のものである。Julia v1.1.1 がリリースされたのが 2019/05/16 であるから,少なくともそれより前のバージョンということになる。

しかしその後の Julia の発展が著しく,今では廃止された関数があったりして,そのままでは参考にならない。また,示された一元配置分散分析は各群の例数が等しく,また,各群の等分散を仮定する検定なので,等分散を仮定しない検定も示しておく。

廃止あるいは仕様変更のあった関数

dat = CSV.read("foo.csv", allowmissing=:none)

dat = CSV.read("lfoo.csv", DataFrame)

と書くようになった。今でも多くのページで単にファイル名を指定すればよいと書いてあるので,小生も最初は困り果てた。 
また,allowmissing は不要になった(当然のものとして組み込まれた)。

by() によるグループごとの統計処理は廃止された。

by(dat, :parasite, :growth => mean)

は groupby() と combine() で以下のように書く。 combine() は非常に強力である。

gd = groupby(dat, :parasite);
combine(gd, :growth => mean)
stat = combine(gd, nrow => :sample_size, :growth_rate => mean => :group_mean, :growth_rate => var => :group_var)

名前の変更は names!() ではなく rename!(),
列の削除は deletecols!() ではなく,select!()
列の並べ替えも permutecols!() ではなく,select!() で行える。

さて,以上の点に留意して,以下のような処理手順を示す。

まずは,データ入力から前処理まで。

using CSV, DataFrames, Statistics, Distributions

dat = CSV.read("Daphniagrowth.csv", DataFrame);
# '.' を使っている列名は,最初に '_' に変更しておく
rename!(dat, "growth.rate" => "growth_rate");
# 群ごとの統計は groupby(), combine() を使う
gd = groupby(dat, :parasite);
stat = combine(gd, nrow => :sample_size, :growth_rate => mean => :group_mean, :growth_rate => var => :group_var)
println(stat)

4×4 DataFrame

 Row │ parasite                   sample_size  group_mean  group_var 
     │ String                     Int64        Float64     Float64   
─────┼───────────────────────────────────────────────────────────────
   1 │ control                             10    1.21391   0.0130376
   2 │ Metschnikowia bicuspidata           10    0.801154  0.0465879
   3 │ Pansporella perplexa                10    1.07636   0.0322315
   4 │ Pasteuria ramosa                    10    0.482203  0.0375737

後は,計算処理を黙々と行う(ドット演算子に注意)
まずは,各群の分散が等しいことを仮定する検定。

k = length(gd)
sample_size = stat[!, :sample_size];
group_mean  = stat[!, :group_mean ];
group_var   = stat[!, :group_var  ];
n = sum(sample_size)
grand_mean = sum(sample_size .* group_mean) / n
St = var(dat[!, :growth_rate]) * (n - 1)
Sb = sum(sample_size .* (group_mean .- grand_mean) .^ 2)
Sw = St - Sb
dfb = k - 1
dfw = n - k
dft = n - 1
Vb = Sb / dfb
Vw = Sw / dfw
# 最終的に必要なのは,F 統計量と,P 値
F = Vb / Vw
p_value = ccdf(FDist(dfb, dfw), F)
# 見やすくするためにデータフレームにしておこう
ANOVA_table = DataFrame(SS=[Sb, Sw, St], df=[dfb, dfw, dft],
      MS=[Vb, Vw, ""], F=[F, "", ""], p_value=[p_value, "", ""])
println(ANOVA_table)

3×5 DataFrame

 Row │ SS       df     MS         F        p_value     
     │ Float64  Int64  Any        Any      Any         
─────┼─────────────────────────────────────────────────
   1 │ 3.13791      3  1.04597    32.3252  2.57128e-10
   2 │ 1.16488     36  0.0323577
   3 │ 4.30278     39


次は,各群の分散が等しいとは仮定しない,いわゆるウェルチの方法によるものである。

grand_mean = mean(dat[!, :growth_rate])
k = length(gd)
sample_size = stat[!, :sample_size];
group_mean = stat[!, :group_mean];
group_var = stat[!, :group_var];
wi = sample_size ./ group_var;
w = sum(wi)
m = sum(wi .* group_mean) / w
a = sum(((1 .- wi ./ w) .^ 2) ./ (sample_size .- 1)) / (k^2 - 1)
F = sum(wi .* (group_mean .- m) .^ 2) / (k - 1) / (1 + 2a*(k - 2))
dfb = k - 1
dfw = 1 / 3a
p_value = ccdf(FDist(dfb, dfw), F)
println("F = $F, p_value = $p_value")

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia の Plots でポリゴンを描画する

2021年01月15日 | ブログラミング

Julia の Plots でポリゴンを描画する方法を探すのに苦労したので,メモしておく。

using Plots

plot(x, y, seriestype=:shape, grid=false, showaxis=false, ticks=false, aspect_ratio=1, label="")

指定が長いので,シュガーコート関数を定義しておくとよいかも。

polygon(x, y) = plot(x, y, seriestype=:shape, grid=false, showaxis=false, ticks=false, aspect_ratio=1, label="")

このようにしておいて,

polygon([1,6,5,3,2], [1,2,6,3,5])

なお,tick ラベルを描かないように ticks=false とするのも,検索には引っかからず,推理してたどり着いた。

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

指定された平均値と標準偏差を持つ正規乱数ベクトルを返す

2021年01月13日 | ブログラミング

平均値 mu,標準偏差 sigma を与え,サンプルサイズ n の 1 変数データを返す

function gendat(n::Int64=1000; mu::Float64=0.0, sigma::Float64=1.0)
        x = randn(n)
        (x .- mean(x)) ./ std(x) .* sigma .+ mu
end

x = gendat(10);
using Statistics
println("mean = $(mean(x)) sd = $(std(x))")

mean = -1.3322676295501878e-16 sd = 1.0

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

julia の関数定義

2021年01月13日 | ブログラミング

引数が違う同じ名前の関数が定義できるのは良いなあ。

function a(x::Vector{Float64})
        return 1
end

function a(n::Int64)
        return 2
end

function a(n::Int64, r::Float64)
        return 3
end


a([1.2, 2.4]) # 1 を返す
a(1)          # 2  を返す
a(1, 0.1)     # 3 を返す

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ガースーさん

2021年01月11日 | 雑感

今すぐ,ですね

緊急事態宣言が必要としている府県の要請に応じてください

外国からの帰国も外国人の入国についても,入国拒否をしてください

できれば,総理大臣を辞任してください

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ネットニュースでも取り上げるのは控えたらよいと思うのだけど

2021年01月11日 | 雑感

ホリエモンとか,桝添さんとか,その他も なんかコンプライアンス的にも不適切な人じゃないの?

こういうことを言ってはいけないのだろうけど,ホリエモンはいわゆる前科者でしょ(^_^;)

桝添さんも怪しいし。

その他には,コンプライアンス的には問題ないかも知れないけど,ひろゆきさんとか,高須さんとか,アッコさんとか,もろもろのいわゆる「ご意見番としての芸能人」。

別に,彼らは誰もが認めているオピニオンリーダーとはいえないと思うんだけど(好きな人はいいんだけど)。

google news では特定の発信元の記事はブロックできるのだけど,livedoor とか auone とかは個別にブロックできないんだよね。

また,数々の発信元をブロックしたあげくだけど,tv 関係はブロックしようとしても,ブロックのメニューさえ出ないんだよね。忖度か...

まあ,ネットニュースなんか見なければいいのだけどね。チャンチャン...

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia で「ハーディ・ラマヌジャン数」を

2021年01月11日 | ブログラミング

昔の記事を自分から見に行くということはほとんどないが,ときどき誰かが読んでくれたという記録に目がとまり,一体どんなことを書いたのかなと確認することがある。

ハーディ・ラマヌジャン数」もそのようなものだった。

それは,R で書いたプログラムであったが,今回 Julia で書いて,速度を比較してみた。

Julia のプログラムは

function ramanujan4()
    m = floor(Int, 9999^(1/3))
    x = (1:m) .^ 3
    z = freqtable([x[i] + x[j] for i=1:m, j=1:m if i < j])
    z[z .>= 2]
end

Benchmark によれば,実行時間の中央値は 56 μs であった。R は 29 ms だったので,Julia は R より,500 倍速いということだ。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia での統計関数についてもまとめておく

2021年01月10日 | ブログラミング

cdf() に対して ccdf(),quantile() に対して cquantile() があるので,修正した。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村