裏 RjpWiki

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

Julia で covid19

2021年03月22日 | ブログラミング

中澤先生の「鐵人三国誌」【第685回】https://minato.sip21c.org/im3r/20210319.htmlでも紹介されている
都道府県別新規感染報告数の推移を描くコード」を Julia でやってみることはできるだろうかと,あれこれやってみた結果。

基本骨格は以下のようになるだろう。

その前に。以下に留意すべし。

We have invested a lot of time and effort in creating COVID-19 Data Hub, please cite the following when using it:

  Guidotti, E., Ardia, D., (2020), "COVID-19 Data Hub", Journal of Open Source Software
  5(51):2376, doi: 10.21105/joss.02376.

A BibTeX entry for LaTeX users is

  @Article{,
    title = {COVID-19 Data Hub},
    year = {2020},
    doi = {10.21105/joss.02376},
    author = {Emanuele Guidotti and David Ardia},
    journal = {Journal of Open Source Software},
    volume = {5},
    number = {51},
    pages = {2376},
  }

さて。では,順をおって記述しましょう。

図を描くので以下を入力。

using  Plots

zip ファイルをダウンロードする関数(country はその後の処理では参照されないようだ)

function doitonlyonece(country = "JPN", level = 2, start = "2010-01-01", raw = true)
  url = "https://storage.covid19datahub.io"
  name = @sprintf("%sdata-%s", (raw ? "raw" : ""), level)
  zip = @sprintf("%s/%s.zip", url, name)
  file = @sprintf("%s.zip", name)
  download(zip, file)
end

最新のデータをダウンロードするなど,
必要なときに 1 回だけ実行する。

zipfile = doitonlyonece()

ダウンロードしたファイルから必要なデータフレームを抽出する関数

using ZipFile, CSV, DataFrames
function readintodataframe(zipfile) # 引数は doitonlyonece() の戻り値
    z = ZipFile.Reader(zipfile)
    csvfile = replace(zipfile, ".zip" => ".csv")
    fileinzip = filter(x->x.name == csvfile, z.files)[1]
    df = CSV.File(read(fileinzip)) |> DataFrame
    df
end

必要なら何回でも実行してかまわない。

df = readintodataframe(zipfile)

描画する都道府県のデータを抽出する(好きなだけどうぞ)。

Tokyo = df[df[:, "administrative_area_level_2"] .== "Tokyo", "confirmed"]
Osaka = df[df[:, "administrative_area_level_2"] .== "Osaka", "confirmed"]
Hyogo = df[df[:, "administrative_area_level_2"] .== "Hyogo", "confirmed"]

図を描く。初回は plot(),二回目以降は plot!()

plot(vcat(Tokyo[1], diff(Tokyo)), label="")
plot!(vcat(Osaka[1], diff(Osaka)), label="")
plot!(vcat(Hyogo[1], diff(Hyogo)), label="")

最後に総括。

plot!()

これらを,もっと使いやすくして提供するのがプログラマーの役目。

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

「都道府県別新規感染報告数の推移を描くコード」に触ってみる

2021年03月22日 | ブログラミング

中澤先生の「鐵人三国誌」【第685回】https://minato.sip21c.org/im3r/20210319.htmlでも紹介されている
都道府県別新規感染報告数の推移を描くコード」を使ってみようと思い,ついでに色々書き足してみたりした。
描画対象とする都道府県を都道府県コードのベクトルで指定したり,片対数グラフの他に普通メモリでも描けるようにしたりが大きな変更である。そのほか,プログラムを関数にしてデータとプログラムをなるべく分離したので,保守しやすくなったと思うがどうだろうか。
呼出方は,以下のようである。描画対象やパラメータを変えて描画するときに,1行だけ入力すればよいので便利かな?

draw.covid19(c(13, 27, 28, 4), event.data)
draw.covid19(11:14, event.data, log="")

 

# Plotting new confirmed cases for Japan against date in semilog scale.
# November 24, 2020: Made
# December 22, 2020: Change color palette
# January 8, 2021: 3 prefs were added
# January 17, 2021: New action date has been added
# February 8, 2021: Emergency is off at Tochigi pref
# February 28, 2021: Emergency is off at several prefs
# March 23, 2021: rewrite easy to use. by 'r-de-r'

draw.covid19 <- function(prefs, event.data,
                         log="y", # plot()
                         bty="n", # legend()
                         lty="dotted", lwd=0.5, col="gray", # segments()
                         lwd2=NULL, lty2=NULL, # lines()
                         pos=4 # text()
                         ) {

    if (!require(COVID19)) {
        install.packages("COVID19", dep=TRUE)
        library(COVID19)
    }

    # 都道府県名(都道府県コード順)
    dict.en <-  c("Hokkaido", "Aomori", "Iwate", "Miyagi", "Akita",
                  "Yamagata", "Fukushima", "Ibaraki", "Tochigi", "Gunma",
                  "Saitama", "Chiba", "Tokyo", "Kanagawa", "Niigata",
                  "Toyama", "Ishikawa", "Fukui", "Yamanashi", "Nagano",
                  "Gifu", "Shizuoka", "Aichi", "Mie", "Shiga",
                  "Kyoto", "Osaka", "Hyogo", "Nara", "Wakayama",
                  "Tottori", "Shimane", "Okayama", "Hiroshima", "Yamaguchi",
                  "Tokushima", "Kagawa", "Ehime", "Kochi", "Fukuoka",
                  "Saga", "Nagasaki", "Kumamoto", "Oita", "Miyazaki",
                  "Kagoshima", "Okinawa")
    dict.jpn <- c("北海道", "青森", "岩手", "宮城",   "秋田",
                  "山形",   "福島", "茨城", "栃木",   "群馬",
                  "埼玉",   "千葉", "東京", "神奈川", "新潟",
                  "富山",   "石川", "福井", "山梨",   "長野",
                  "岐阜",   "静岡", "愛知", "三重",   "滋賀",
                  "京都",   "大阪", "兵庫", "奈良",   "和歌山",
                  "鳥取",   "島根", "岡山", "広島",   "山口",
                  "徳島",   "香川", "愛媛", "高知",   "福岡",
                  "佐賀",   "長崎", "熊本", "大分",   "宮崎",
                  "鹿児島", "沖縄")

    # ラベル,タイトル
    xlab <- "時間経過(年)"
    ylab <- "新規確定患者報告数"
    main <- "COVID-19の都道府県別新規確定患者報告数の推移
               [Ref.] Guidotti E, Ardia D (2020) COVID-19 Data Hub.
              Working paper https://doi.org/10.13140/RG.2.2.11649.81763"

    # 描画範囲(横軸)
    from <- sprintf("2020-%02d-01", 1:12)
    to   <- sprintf("2021-%02d-01", 1:3)

    # 発生事象描画関数
    event.draw <- function(x) {
        factor <- function(x) if (log == "") log(x)*350 - 250 else x
        date.type = as.Date(x[1,1])
        segments(date.type, 1, date.type, factor(8000), lty=lty, lwd=lwd, col=col)
        text(date.type, factor(x[1,3]), x[1,2],
            cex=0.8, offset=0, pos=pos, xpd=TRUE)
    }

    options(scipen=7)

    # RStudio/Mac OS の場合に文字化け回避のために必要
    par(family= "Hiragino Mincho Pro W3") # ほかに "HiraKakuProN-W3", "Hiragino Kaku Gothic Pro W3"

    # 折線の太さ・線種の定義
    np <- length(prefs)
    if (is.null(lwd2) || length(lwd2) != np) {
        lwd2 <- rep(0.5, np)
    }
    if (is.null(lty2) || length(lty2) != np) {
        lty2 <- rep("solid", np)
    }

    # 都道府県名取り出し
    prefsJ <- dict.jpn[prefs]
    prefs  <- dict.en [prefs]

    # cols <- rainbow(length(prefs))
    cols <- 1:length(prefs)
    palette("Okabe-Ito")
    y <- covid19(country="JPN", level=2, start="2020-01-01", end=Sys.Date())
    y <- subset(y, date < max(date))
    y$prefs <- as.factor(y$administrative_area_level_2)
    cfd <- c(y$confirmed[1], diff(y$confirmed))
    cfd <- ifelse(cfd < 1, NA, cfd)
    plot(cfd ~ as.Date(date, "%d/%m/%y"), data=y,
        xlim=c(min(y$date), max(y$date) + 20), ylim=c(1, ifelse(log=="", 2800, 10000)),
        axes=FALSE, log=log, type="n", frame=FALSE,
        xlab=xlab, ylab=ylab, main=main,
        sub=sprintf("%s まで", max(y$date)))
    axis(1, c(as.Date(from), as.Date(to)), c(from, to))
    if (log == "y") {
        axis(2, 10^(0:4))
    } else {
        axis(2)
    }
    for (i in 1:nrow(event.data)) {
        event.draw(event.data[i,])
    }
    for (i in levels(y$prefs)) {
        dd <- subset(y, prefs == i)
        dd$confirmedpp <- c(dd$confirmed[1], diff(dd$confirmed))
        dd$confirmedpp <- ifelse(dd$confirmedpp < 0, 0, dd$confirmedpp)
        if (i %in% prefs) {
            index <- grep(i, prefs)
            lines(dd$date, dd$confirmedpp, col=cols[index],
                lwd=lwd2[index], lty=lty2[index])
            text(dd$date[length(dd$date)]+1,
                max(dd$confirmedpp[length(dd$confirmedpp)], 1),
                prefsJ[index], col=cols[index], pos=4, xpd=TRUE)
        }
    }
    legend("topleft", col=cols, lwd=2, legend=prefsJ, bty=bty)
}

# 発生事象データ(随時追加)
event.data = data.frame(
    date =  c("2020-04-08", "2020-05-14", "2020-06-19", "2020-07-22", "2020-10-01",
              "2021-01-07", "2021-01-14", "2021-02-08", "2021-02-28"),
    event = c("緊急事態宣言発出", "39県解除", "東京解除", "GoToキャンペーン",
              "GoTo東京追加,\nGoToイート開始", "一都三県\n緊急事態宣言発出",
              "二府五県\n緊急事態宣言追加", "栃木県\n緊急事態宣言解除",
              "二府四県\n緊急事態宣言解除"),
    pos =   c(5000, 3000, 1700, 5000, 3000,
              5000, 3000, 1700, 1000))

# 描画都道府県は draw.covid19() の第1パラメータで c(13, 27, 4) のように都道府県コードで指定する
# 折線グラフの詳細指定(省略可)
lwd2 = rep(0.5, 4); lwd2[c(1, 4)] = 2 # 太さ
lty2 = rep("solid", 4)          # 線種

draw.covid19(c(13, 27, 28, 4), event.data, lwd2=lwd2, lty2=lty2)
draw.covid19(c(13, 27, 28, 4), event.data, log="", lwd2=lwd2, lty2=lty2)

 

 

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

Julia に翻訳--116 サーストンの一対比較法

2021年03月22日 | ブログラミング

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

サーストンの一対比較法
http://aoki2.si.gunma-u.ac.jp/R/ThurstonePairedComparison.html

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

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

関数名が長い。

==========#

using Rmath, Statistics, Printf, Plots

function thurstonepairedcomparison(x; digits=5)
    nr, nc = size(x)
    nr == nc || error("required square matrix")
    n = [i == j ? 1 : x[i, j] + x[j, i] for i = 1:nr, j = 1:nr]
    x = [i == j ? missing : qnorm(x[i, j] / n[i, j]) for i = 1:nc, j = 1:nc]
    score = [mean(skipmissing(x[i, :])) for i = 1:nc]
    sortedscore = sort(score)
    printthurstonepairedcomparison(score, sortedscore, digits)
    plotthurstonepairedcomparison(score, nc)
    Dict(:score => score, :sortedscore => sort(score))
end

function printthurstonepairedcomparison(score, sortedscore, digits)
    println("スコア")
    println(round.(score, digits=digits))
    println("\nソートされたスコア")
    println(round.(sortedscore, digits=digits))
end

function plotthurstonepairedcomparison(score, nc)
    pyplot(size=(600, 60))
    plt = scatter(score, repeat([25], nc), grid=false, tick_direction=:out,
            yshowaxis=false, ylims=(0, 60), yticks=false, label="")
    for i = 1:nc
        annotate!(score[i], 50, text(string(i), 10, :red), label="")
    end
    display(plt)
end

x = [ 0 72 73 70 73 73 69 60
      2  0 32  1 18  6  2  1
      1 42  0  1 19 17  2  1
      4 73 73  0 59 63 44 22
      1 56 55 15  0 24  9  7
      1 68 57 11 50  0 15 10
      5 72 72 30 65 59  0 18
     14 73 73 52 67 64 56  0 ];

thurstonepairedcomparison(x)
# スコア
# [1.79165, -1.50557, -1.39753, 0.62813, -0.66134, -0.37066, 0.4889, 1.02641]

# ソートされたスコア
# [-1.50557, -1.39753, -0.66134, -0.37066, 0.4889, 0.62813, 1.02641, 1.79165]

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

Julia に翻訳--115 対応のあるデータの二つの相関係数の相等性の検定

2021年03月22日 | ブログラミング

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

対応のあるデータの二つの相関係数の相等性の検定
http://aoki2.si.gunma-u.ac.jp/R/diff_r2.html

ファイル名: diffr2.jl
関数名:    diffr2(データ行列)
          diffr2(相関係数行列, サンプルサイズ)

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

データ行列を与える場合と,cor と n を与える場合の 2 通りの関数を定義した。

==========#

using Statistics, Rmath

function diffr2(x)
    n, nc = size(x)
    nc == 4 && n != nc || error("required nx4 data matrix")
    diffr2(cor(x), n)
end

function diffr2(r, n)
    z12 = atanh(r[1,2])
    z34 = atanh(r[3,4])
    a1 = r[1, 3] * r[2, 4] + r[1, 4] * r[2, 3]
    a2 = -r[3, 4] * (r[1, 3] * r[2, 3] + r[1, 4] * r[2, 4])
    a3 = -r[1, 2] * (r[1, 3] * r[1, 4] + r[2, 3] * r[2, 4])
    a4 = r[1, 2] * r[3, 4] * (r[1, 3]^2 + r[1, 4]^2 + r[2, 3]^2 + r[2, 4]^2) / 2
    d = (1 - r[1, 2]^2) * (1 - r[3, 4]^2)
    chisq = (n - 3) * (z12 - z34)^2 / (2 - 2 * (a1 + a2 + a3 + a4) / d)
    df = 1
    p = pchisq(chisq, df, false)
    println("chisq. = $chisq,  df = $df,  p value = $p")
    Dict(:chisq => chisq, :df => df, :pvalue => p)
end

include("gendat.jl")
d = gendat(220, [0.439,0.288,0.329,0.354,0.320,0.595]);
cor(d)
# 4×4 Array{Float64,2}:
#  1.0    0.439  0.288  0.329
#  0.439  1.0    0.354  0.32
#  0.288  0.354  1.0    0.595
#  0.329  0.32   0.595  1.0
diffr2(d)
# chisq. = 5.500229788108368,  df = df,  p value = 0.01901397495593562

r = [  1.0    0.439  0.288  0.329
       0.439  1.0    0.354  0.32
       0.288  0.354  1.0    0.595
       0.329  0.32   0.595  1.0   ];
diffr2(r, 220)
# chisq. = 5.500229788108428,  df = 1,  p value = 0.019013974955934914

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

Julia に翻訳--114 同じサンプルからの相関係数の差

2021年03月22日 | ブログラミング

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

同じサンプルからの相関係数の差
http://aoki2.si.gunma-u.ac.jp/R/diff_r.html

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

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

Julia では t という変数名が許されるので気持ちがよい。

==========#

using Rmath

function diffr(n, rxy, rvy, rxv)
    detR = (1 - rxy^2 - rvy^2 - rxv^2) + 2 * rxy * rxv * rvy
    t = (abs(rxy - rvy) * sqrt((n - 1) * (1 + rxv))) /
        sqrt(2 * detR * (n - 1) / (n - 3) + (rxy + rvy)^2 * (1 - rxv)^3 / 4)
    df = n - 3
    p = pt(t, df, false) * 2
    println("t = $t,  df = $df,  p value = $p")
    Dict(:t => t, :df => df, :pvalue => p)
end

diffr(50, 0.5, 0.32, 0.65)
# t = 1.696399559067307,  df = 47,  p value = 0.09642490775716749

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

Julia に翻訳--113 標本相関係数の同等性の検定

2021年03月22日 | ブログラミング

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

標本相関係数の同等性の検定
http://aoki2.si.gunma-u.ac.jp/R/eq-cor.html

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

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

==========#

using Rmath

function eqcor(n, r)
    (all(n .> 3) && length(n) == length(r) && all(floor.(n) .== n) && all(abs.(r) .<= 1)) ||
        error("parameter error")
    k = length(n)
    v = n .- 3
    z = atanh.(r)
    sv = sum(v)
    svz = sum(v .* z)
    chi = sum(v .* z .^ 2) - svz^2 / sv
    df = k - 1
    pvalue = pchisq(chi, df, false)
    println("chisq. = $chi,  df = $df,  p value = $pvalue")
    dic = Dict(:chisq => chi, :df => df, :pvalue => pvalue)
    if pvalue > 0.05
        est = tanh(svz / sv)
        println("estimated rho = $est")
        dic[:est] = est
    end
    dic
end

n = [10, 16, 8, 29, 36]
r = [0.658, 0.285, 0.569, 0.427, 0.374]
eqcor(n, r)
# chisq. = 1.4239560404291218,  df = 4,  p value = 0.8400208372208527
# estimated rho = 0.4179631144001677

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

Julia に翻訳--112 母相関係数が 0 以外の特定の値であるかどうかの検定

2021年03月22日 | ブログラミング

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

母相関係数が 0 以外の特定の値であるかどうかの検定
http://aoki2.si.gunma-u.ac.jp/R/corr.html

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

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

短!

==========#

using Rmath

function rtest(n, r, rho)
    z = abs(atanh(r) - atanh(rho)) * sqrt(n - 3)
    p = pnorm(z, false) * 2
    println("z = $z,  p value = $p")
    Dict(:z => z, :pvalue => p)
end

rtest(24, 0.476, 0.3)
# z = 0.954459123362117,  p value = 0.33985129105621764

 

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

Julia に翻訳--111 相関係数の有意性検定

2021年03月22日 | ブログラミング

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

相関係数の有意性検定
http://aoki2.si.gunma-u.ac.jp/R/cor2.html

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

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

==========#

using Rmath

function cor2test(n, r; conflevel=0.95, method = "pearson") # "kendall", "spearman"
    if method != "kendall"
        t = abs(r) * sqrt((n - 2) / (1 - r^2))
        df = n - 2
        p = pt(t, df, false) * 2
        q = qnorm(0.5 - conflevel / 2) / sqrt(n - 3)
        confint = tanh.(atanh(r) .+ [q, -q])
        println("t = $t,  df = $df,  p = $p\nconfint = $confint")
        Dict(:t => t, :df => df, :pvalue => p, :confint => confint)
    else
        z = abs(r) / sqrt((4 * n + 10) / (9 * n * (n - 1)))
        p = pnorm(z, false) * 2
        println("z = $z,  p = $p")
        Dict(:Z => z, :pvalue => p)
    end
end


cor2test(24, 0.476)
# t = 2.538688822593072,  df = 22,  p = 0.018712499984981517
# confint = [0.08985742931333052, 0.7377384989803106]

cor2test(24, 0.282, method="spearman")
# t = 1.378650599523675,  df = 22,  p = 0.18186096371901114
# confint = [-0.13697916849860975, 0.6153911248557894]

cor2test(24, 0.493, method="kendall")
# z = 3.3750855083507507,  p = 0.0007379275879158393

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

Julia に翻訳--110 共分散分析

2021年03月22日 | ブログラミング

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

共分散分析
http://aoki2.si.gunma-u.ac.jp/R/covar-test.html

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

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

==========#

using Statistics, Rmath, Printf

function covartest(x, y, g)
    index, nj = table(g)
    n = sum(nj)
    k = length(nj)
    mx = mean(x)
    mxj = [mean(x[g .== g0]) for g0 in index]
    sstx = (n - 1) * var(x)
    ssbx = sum(nj .* (mxj .- mx) .^ 2)
    sswx = sstx - ssbx
    my = mean(y)
    myj = [mean(y[g .== g0]) for g0 in index]
    ssty = (n - 1) * var(y)
    ssby = sum(nj .* (myj .- my) .^ 2)
    sswy = ssty - ssby
    spt = (n - 1) * cov(x, y)
    spb = sum(nj .* (mxj .- mx) .* (myj .- my))
    spw = spt - spb
    sswy = sswy - spw^2 / sswx
    ssty = ssty - spt^2 / sstx
    ssby = ssty - sswy
    hensax = x .- mxj[g]
    hensay = y .- myj[g]
    xy = hensax .* hensay
    xx = hensax .^ 2
    numerator = [sum(xy[g .== g0]) for g0 in index] # tapply(xy, g, sum)
    denominator = [sum(xx[g .== g0]) for g0 in index] # tapply(xx, g, sum)
    a = numerator ./ denominator
    b = myj .- a .* mxj
    predicty = a[g] .* x .+ b[g]
    sswyj = sum((y .- predicty) .^ 2)
    dfr = k - 1
    dfe = n - 2 * k
    dft = n - k - 1
    diffreg = sswy - sswyj
    msr = diffreg / dfr
    mse = sswyj / dfe
    msw = sswy / dft
    f = msr / mse
    p = pf(f, dfr, dfe, false)
    part1="H0: 各群の回帰直線の傾きは同じである"
    names = ["group x slope", "residual", "total"]
    ss = [diffreg, sswyj, sswy]
    df = [dfr, dfe, dft]
    ms = [msr, mse, msw]
    anovatable(part1, names, ss, df, ms, f, dfr, dfe, p)
    if p > 0.05
        msby = ssby / dfr
        mswy = sswy / dft
        f2 = msby / mswy
        p2 = pf(f2, dfr, dft, false)
        part2 = "H0: 共変量で調整した平均値は同じである"
        names2 = ["effect & group", "residual", "total"]
        ss2 = [ssby, sswy, ssty]
        df2 = [dfr, dft, n-2]
        ms2 = [msby, mswy, ssty/(n-2)]
        println()
        anovatable(part2, names, ss2, df2, ms2, f2, dfr, dft, p2)
    end
end

function table(x) # indices が少ないとき
    indices = sort(unique(x))
    counts = zeros(Int, length(indices))
    for i in indexin(x, indices)
        counts[i] += 1
    end
    return indices, counts
end

function anovatable(title, names, ss, df, ms, f, df1, df2, p)
    println(title)
    @printf("%14s  %12s  %3s  %12s\n", "", "SS", "df", "MS")
    for i = 1:3
        @printf("%14s  %12.7f  %3d  %12.7f\n", names[i], ss[i], df[i], ms[i])
    end
    println("F value = $f,  df1 = $df1,  df2 = $df2,  p value = $p")
end

dat = [ 5.9 1.87 1
        5.7 1.19 1
        5.7 1.22 1
        6.7 1.46 1
        6.5 1.35 1
        6.2 1.16 1
        6.3 1.62 1
        6.7 2.00 1
        5.6 1.28 2
        5.3 1.06 2
        5.5 1.17 2
        5.6 1.74 2
        5.5 1.13 2
        5.3 1.18 2
        5.2 1.03 2
        5.6 1.23 2
        5.5 0.90 2
        5.5 1.24 2];
x = dat[:, 1];
y = dat[:, 2];
g = Int.(dat[:, 3]); # 整数型であること!
covartest(x, y, g)

# H0: 各群の回帰直線の傾きは同じである
#                           SS   df            MS
#  group x slope     0.0357661    1     0.0357661
#       residual     0.9187879   14     0.0656277
#          total     0.9545539   15     0.0636369
# F value = 0.544984285639864,  df1 = 1,  df2 = 14,  p value = 0.472568833125542

# H0: 共変量で調整した平均値は同じである
#                           SS   df            MS
#  group x slope     0.0000079    1     0.0000079
#       residual     0.9545539   15     0.0636369
#          total     0.9545618   16     0.0596601
# F value = 0.00012428977991377308,  df1 = 1,  df2 = 15,  p value = 0.9912518692175619

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

Julia に翻訳--109 多次元分布の平均値の差の検定,ウィルクスのΛ

2021年03月22日 | ブログラミング

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

多次元分布の平均値の差の検定(ウィルクスのΛ)
http://aoki2.si.gunma-u.ac.jp/R/wilks.html

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

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

reduce(+, vars, dims=3) の結果も 3 次元(n x n x 1 だけど)とは?
ちょっと気持ち悪いので,単純に結果を累積するようにプログラムを書き換える。

==========#

using Statistics, Rmath, LinearAlgebra

function wilks(dat, g)
    nv = size(dat, 2)
    index, gcase = table(g)
    n = sum(gcase)
    k = length(gcase)
    w = zeros(nv, nv)
    for (i, gvar) in enumerate(index)
        z = dat[g .== gvar, :]
        w += cov(z) * (size(z, 1) - 1)
    end
    t = cov(dat) * (n - 1)
    Λ = det(w) / det(t)
    sl = sqrt(Λ)
    if nv == 2
        f = (1 - sl) * (n - k - 1) / sl / (k - 1)
        df1 = 2 * (k - 1)
        df2 = 2 * (n - k - 1)
        p = pf(f, df1, df2, false)
        println("F = $f,  df1 = $df1,  df2 = $df2,  p vallue = $p")
        Dict(:Fvalue => f, :df1 => df1, :df2 => df2, :pvalue => p)
    elseif k == 2
        f = (1 - Λ) * (n - k - nv + 1) / Λ / nv
        df2 = n - k - nv + 1
        p = pf(f, nv, df2, false)
        println("F = $f,  df1 = $nv,  df2 = $df2,  p vallue = $p")
        Dict(:Fvalue => f, :df1 => nv, :df2 => df2, :pvalue => p)
    elseif k == 3
        f = (1 - sl) * (n - k - nv + 1) / sl / nv
        df1 = 2 * nv
        df2 = 2 * (n - k - nv + 1)
        p = pf(f, df1, df2, false)
        println("F = $f,  df1 = $df1,  df2 = $df2,  p vallue = $p")
        Dict(:Fvalue => f, :df1 => df1, :df2 => df2, :pvalue => p)
    else
        chisq = ((nv + k) / 2 - n + 1) * log(Λ)
        df = nv * (k - 1)
        p = pchisq(chisq, df, false)
        println("chisq. = $chisq,  df = $df,  p vallue = $p")
        Dict(:chisq => chisq, :df => df, :pvalue => p)
    end
end

function table(x) # indices が少ないとき
    indices = sort(unique(x))
    counts = zeros(Int, length(indices))
    for i in indexin(x, indices)
        counts[i] += 1
    end
    return indices, counts
end

dat0 = [1 2.9 161.7 120.8
        1 2.3 114.8 85.2
        1 2 128.4 92
        1 3.2 149.2 97.3
        1 2.7 126 81.1
        1 4.4 133.8 107.6
        1 4.1 161.3 114
        1 2.1 111.5 77.3
        2 4.8 198.7 172.9
        2 3.6 199.3 157.9
        2 2 188.4 152.7
        2 4.9 183.6 164.2
        2 3.9 173.5 172.2
        2 4.4 184.9 163.2];
dat = dat0[:, 2:4];
g = dat0[:, 1];
wilks(dat, g)
# F = 33.80499798162517,  df1 = 3,  df2 = 10,  p vallue = 1.516644831663382e-5

wilks(dat0[:, 2:3], dat0[:, 1])
# F = 10.911143164354023,  df1 = 2,  df2 = 22,  p vallue = 0.0005105098610968934

dat2 = vcat(dat0, [3 2.3 195.2 154.4; 3 3.6 189.2 150.8; 3 3.4 179.3 182.4])
wilks(dat2[:, 2:4], dat2[:, 1])
# F = 9.531929948587386,  df1 = 6,  df2 = 24,  p vallue = 2.143115741229085e-5

dat3 = vcat(dat2, [4 2.5 165.2 157.1; 4 4.6 192.2 149.8; 4 5.2 181.2 180.6])
wilks(dat3[:, 2:4], dat3[:, 1])
# chisq. = 38.2334484265828,  df = 9,  p vallue = 1.5827976779857345e-5

 

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

Julia に翻訳--108 三要因の分散分析,SABC タイプ,RBFpqr デザイン,被検者内計画

2021年03月22日 | ブログラミング

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

三要因の分散分析(SABC タイプ;RBFpqr デザイン;被検者内計画)
http://aoki2.si.gunma-u.ac.jp/R/SABC.html

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

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

sum(tbl, dims=) のオンパレード

==========#

using Rmath, Statistics, Printf

function sabc(tbl)
    r, q, p, n = size(tbl)
    X = sum(tbl)^2/n/p/q/r
    A = sum(sum(tbl, dims=[1,2,4]) .^ 2) / n / q / r
    B = sum(sum(tbl, dims=[1,3,4]) .^ 2) / n / p / r
    C = sum(sum(tbl, dims=[2,3,4]) .^ 2) / n / p / q
    S = sum(sum(tbl, dims=[1,2,3]) .^ 2) / p / q / r
    AB = sum(sum(tbl, dims=[1,4]) .^ 2) / n / r
    AC = sum(sum(tbl, dims=[2,4]) .^2) / n / q
    BC = sum(sum(tbl, dims=[3,4]) .^2) / n / p
    AS = sum(sum(tbl, dims=[1,2]) .^ 2) / q / r
    BS = sum(sum(tbl, dims=[1,3]) .^ 2) / p / r
    CS = sum(sum(tbl, dims=[2,3]) .^2) / p / q
    ABS = sum(sum(tbl, dims=1) .^ 2) / r
    ACS = sum(sum(tbl, dims=2) .^ 2) / q
    BCS = sum(sum(tbl, dims=3) .^ 2) / p
    ABC = sum(sum(tbl, dims=4) .^ 2) / n
    ABCS = sum(tbl .^ 2)
    SSA = A - X
    SSB = B - X
    SSC = C - X
    SSAB = AB - A - B + X
    SSAC = AC - A - C + X
    SSBC = BC - B - C + X
    SSABC = ABC - AB - AC - BC + A + B + C - X
    SST = ABCS - X
    SSS = S - X
    SSAS = AS - A - S + X
    SSBS = BS - B - S + X
    SSCS = CS - C - S + X
    SSABS = ABS - AB - AS - BS + A + B + S - X
    SSACS = ACS - AC - AS - CS + A + C + S - X
    SSBCS = BCS - BC - BS - CS + B + C + S - X
    SSABCS = ABCS - ABC - ABS - ACS - BCS +
             AB + AC + BC + AS + BS + CS - A - B - C - S + X
    SS = [SSS, SSA, SSAS, SSB, SSBS, SSC, SSCS, SSAB,
          SSABS, SSAC, SSACS, SSBC, SSBCS, SSABC, SSABCS, SST]
    n1 = n - 1
    p1 = p - 1
    q1 = q - 1
    r1 = r - 1
    df = [n1, p1, p1 * n1, q1, q1 * n1, r1, r1 * n1, p1 * q1,
          p1 * q1 * n1, p1 * r1, p1 * r1 * n1, q1 * r1, q1 * r1 * n1,
          p1 * q1 * r1, p1 * q1 * r1 * n1, n * p * q * r - 1]
    MS = SS ./ df
    suf = 2:2:14
    Fvalue = fill(NaN, 16);
    pvalue = fill(NaN, 16);
    Fvalue[suf] = MS[suf] ./ MS[suf .+ 1];
    pvalue[suf] = pf.(Fvalue[suf], df[suf], df[suf .+ 1], false)
    rownames = anovatable(SS, df, MS, Fvalue, pvalue)
    names, means, SE, LCL, UCL = summarysabc(r, q, p, n, tbl)
    Dict(:rownames => rownames, :SS => SS, :df => df, :MS => MS,
         :Fvalue => Fvalue, :pvalue => pvalue,
         :names => names, :means => means, :SE => SE, :LCL => LCL, :UCL => UCL)
end

function anovatable(SS, df, MS, Fvalue, pvalue)
    rownames = ["S", "A", "AxS", "B", "BxS", "C", "CxS", "AxB", "AxBxS",
                "AxC", "AxCxS", "BxC", "BxCxS", "AxBxC", "AxBxCxS", "T"]
    println("分散分析表")
    @printf("%7s  %11s  %3s  %12s  %12s  %9s\n",
            "", "SS", "df", "MS", "F value", "p value")
    for i = 1:16
        if i <= 14 && i % 2 == 0
            @printf("%7s  %11.7f  %3d  %12.7f  %12.7f  %9.5f\n",
                rownames[i], SS[i], df[i], MS[i], Fvalue[i], pvalue[i])
        else
            @printf("%7s  %11.7f  %3d  %12.7f\n",
                rownames[i], SS[i], df[i], MS[i])
        end
    end
    rownames
end

function summarysabc(r, q, p, n, tbl)
    meansa = reshape(sum(tbl, dims=[1,2]), 2, 4) / q / r;
    meansb = reshape(sum(tbl, dims=[1,3]), 2, 4) / p / r;
    meansc = reshape(sum(tbl, dims=[2,3]), 3, 4) / p / q;
    means = vcat(mean(meansa, dims=2), mean(meansb, dims=2), mean(meansc, dims=2))
    SE = vcat(std(meansa, dims=2), std(meansb, dims=2), std(meansc, dims=2)) / sqrt(n)
    t95 = qt(0.975, n - 1)
    LCL = means .- t95 .* SE
    UCL = means .+ t95 .* SE
    names = vcat(" a=" .* string.(1:p), (" b=" .* string.(1:q)), (" c=" .* string.(1:r)))
    println("\n推定周辺平均")
    @printf("%3s  %10s  %10s  %10s  %10s\n", "", "Mean", "SE", "LCL", "UCL")
    for i = 1:length(names)
        @printf("%3s  %10.6f  %10.6f  %10.6f  %10.6f\n",
            names[i], means[i], SE[i], LCL[i], UCL[i])
    end
    names, means, SE, LCL, UCL
end

dat = [2,5,9, 6,3,6, 1,2,1, 5,5,5,
  6,7,10, 6,6,7, 2,1,5, 3,6,5,
  5,9,13, 8,5,5, 5,4,3, 4,6,9,
  7,9,14, 10,8,6, 2,5,5, 6,7,7];
tbl = reshape(dat, 3, 2, 2, 4);
sabc(tbl)

# 分散分析表
#                   SS   df            MS       F value    p value
#       S   60.3333333    3    20.1111111
#       A   96.3333333    1    96.3333333    51.0000000    0.00565
#     AxS    5.6666667    3     1.8888889
#       B    3.0000000    1     3.0000000     1.4210526    0.31893
#     BxS    6.3333333    3     2.1111111
#       C   33.5000000    2    16.7500000    60.3000000    0.00011
#     CxS    1.6666667    6     0.2777778
#     AxB   56.3333333    1    56.3333333   101.4000000    0.00209
#   AxBxS    1.6666667    3     0.5555556
#     AxC    6.1666667    2     3.0833333     2.9210526    0.13007
#   AxCxS    6.3333333    6     1.0555556
#     BxC   24.5000000    2    12.2500000     5.3780488    0.04591
#   BxCxS   13.6666667    6     2.2777778
#   AxBxC   41.1666667    2    20.5833333     6.0737705    0.03614
# AxBxCxS   20.3333333    6     3.3888889
#       T  377.0000000   47     8.0212766

# 推定周辺平均
#            Mean          SE         LCL         UCL
# a=1    7.166667    0.790569    4.650722    9.682611
# a=2    4.333333    0.540062    2.614616    6.052051
# b=1    5.500000    0.819327    2.892537    8.107463
# b=2    6.000000    0.504608    4.394111    7.605889
# c=1    4.875000    0.616610    2.912671    6.837329
# c=2    5.500000    0.743023    3.135369    7.864631
# c=3    6.875000    0.599479    4.967190    8.782810

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

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

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