裏 RjpWiki

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

どうやって入力する? 〃 その他

2021年04月23日 | 雑感

これは,々とおなじで「おなじ」といれて変換する。

他に「おなじ」は,ゝ,ゞ,ヽ,ヾ,仝 がある。

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

どうやって入力する? 之

2021年04月23日 | 雑感

これは簡単。

「これ」っていれて,変換する。名前によく使われるので「ゆき」のほうがよく使われているかも。「し」でもよいけど,これは場合によっては変換キーを何回も押さなくてはならなくて,行きすぎたりする。

石碑などの刻印に,「○○建之」とあるのは「○○がこれを建てた」ということ。

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

Julia の小ネタ--19 文字列の連結法 3 通り

2021年04月23日 | ブログラミング

文字列の連結 3 通り

greet = "Hello"
whom = "world"

から,"Hello, world.\n" を作る

第1の方法
string(greet, ", ", whom, ".\n")

第2の方法
greet * ", " * whom * ".\n"

第3の方法
"$greet, $whom.\n"

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

Julia に翻訳--199 因子分析の適合度検定

2021年04月23日 | ブログラミング

#==========
Julia の修行をするときに,いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

因子分析の適合度検定
http://aoki2.si.gunma-u.ac.jp/R/fa.fit.test.html

ファイル名: fafittest.jl  関数名: fafittest

翻訳するときに書いたメモ

計算の途中で結果が Matrix{Any} になることがある。今回の場合は

sc = [i == j ? 1 / sqrt(uniquenesses[i]) : 0 for i = 1:p, j = 1:p]

では,sc が Matrix{Real} になるのが原因。それを回避するには

sc = [i == j ? 1.0 / sqrt(uniquenesses[i]) : 0.0 for i = 1:p, j = 1:p]

とすればよいようだが,なんだかなあ!

==========#

using DataFrames, Statistics, StatsBase, Rmath, LinearAlgebra

function fafittest(data, factors, uniquenesses)
    n, p = size(data)
    r = cor(data)
    sc = [i == j ? 1.0 / sqrt(uniquenesses[i]) : 0.0 for i = 1:p, j = 1:p]
    r = sc * r * sc
    values = eigvals(r) # Julia ではデフォルトで昇順
    e = values[1:p - factors]
    s = -sum(log.(e) .- e) + factors - p
    chisq = (n - 1 - (2 * p + 5) / 6 - (2 * factors) / 3) * s
    df = ((p - factors) ^ 2 - p - factors) / 2
    pvalue = pchisq(chisq, df, false)
    println("chisq. = $chisq,  df = $df,  p value = $pvalue")
    Dict(:chisq => chisq, :df => df, :pvalue => pvalue)
end

data = [
  935  955  926  585 1010  925 1028  807  769   767
  817  905  901  632 1004  950  957  844  781   738
  768  825  859  662  893  900  981  759  868   732
  869  915  856  448  867  874  884  802  804   857
  787  878  880  592  871  874  884  781  782   807
  738  848  850  569  814  950  957  700  870   764
  763  862  839  658  887  900 1005  604  709   753
  795  890  841  670  853  874  859  701  680   772
  903  877  919  460  818  849  884  700  718   716
  761  765  881  485  846  900  981  728  781   714
  747  792  800  564  796  849  932  682  746   767
  771  802  840  609  824  874  859  668  704   710
  785  878  805  527  911  680  884  728  709   747
  657  773  820  612  810  849  909  698  746   771
  696  785  791  578  774  725  932  765  706   795
  724  785  870  509  746  849  807  763  724   760
  712  829  838  516  875  725  807  754  762   585
  756  863  815  474  873  725  957  624  655   620
  622  759  786  619  820  769  807  673  698   695
  668  753  751  551  834  849  807  601  655   642];

uniquenesses = [0.005, 0.186849474704386, 0.235609633298087,
  0.594548313805094, 0.005, 0.005, 0.644059344661751,
  0.005, 0.466720965448361, 0.573892982555005]
fafittest(data, 4, uniquenesses)
# chisq. = 15.163572526523414,  df = 11.0,  p value = 0.1751303147464357

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

Julia の小ネタ--018 演算子 \

2021年04月23日 | ブログラミング

どこかでだれかが,「演算子 \ なんて使い道あるのか?」なんて書いていたが,ありますよ。

a*x = b のとき x を求める,つまり x = b /a = 1/a * b という記述を,\ を使って x = a \ b と記述できるということですよ。つまり,a の逆数と b の積を意味するんですね。

julia> a = 3
3

julia> b = 7
7

julia> x = b / a
2.3333333333333335

julia> x = a \ b
2.3333333333333335

これだけだとあまりありがたみがないんだけど,a や b が配列やベクトルであっても同じように使えますよ。

A * x = B のとき,x は A の逆行列と B の行列積ですね。上と同じように A の逆数と B の積で x = A^(-1) * B = inv(A) * B です。\ を使えば A \ B な訳です。


julia> A = [2 3; 5 7]
2×2 Matrix{Int64}:
 2  3
 5  7

julia> B = [11, 13]
2-element Vector{Int64}:
 11
 13

julia> inv(A) * B
2-element Vector{Float64}:
 -38.00000000000006
  29.00000000000004

julia> A \ B
2-element Vector{Float64}:
 -38.00000000000006
  29.00000000000004

で,同じになることがわかります。

じゃあ,わざわざ新しい演算子なんか持ち出すことはないだろう,ドッチでもよいじゃないか。

いやいや,そうじゃないんです。

inv(A) * B は A \ B より,3倍ぐらい遅いんです。

とは言っても,8000x8000 の行列 A と要素数 8000 のベクトル B の演算だと,inv(A) * B は 18 秒,A \ B は 5 秒ぐらいなので,五十歩百歩ではありますが。

R ではそれぞれ solve(A)*B ,solve(A, B) に対応するのだけど,同じ条件で,前者は 550 秒,後者は 125 秒ぐらいなので,Julia は相当に速いということではあります。

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

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

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