適合度の検定

 

キーワード

適合度の検定，名義尺度変数，カイ二乗分布，漸近検定，一様性の検定，統計検定，帰無仮説，対立仮説，一様乱数，p値，Julia，度数分布，検定手順，有意水準，第一種の過誤，type I error，棄却，採択，期待度数，検定統計量，有意確率，プログラム，サンプルサイズ

Mac OS の R のオープニングメッセージ

あれは，いつ頃からだったんだろうか...

R のオープニングメッセージが英語になって

Natural language support but running in an English locale

というのが表示されている。エラーメッセージも英語になっている。

いいのだけど。

でも，長い間慣れ親しんだ日本語のオープニングメッセージが懐かしい。

ならば，ターミナルのコマンドラインで

defaults write org.R-project.R force.LANG ja_JP.UTF-8

とする。オープニングメッセージだけでなくエラーメッセージも日本語になる。

英語に戻すには，ターミナルのコマンドラインで

defaults write org.R-project.R force.LANG en_US.UTF-8

とする。

ちなみに，おフランス語なら

defaults write org.R-project.R force.LANG fr_FR.UTF-8

使用可能なロケールhttps://docs.oracle.com/cd/E26924_01/html/E27144/glset.html

 

### 

R のセッション中で一時的に英語にしたり日本語にしたりが必要ならば，

Sys.setenv(LANG = "en_US.UTF-8")
Sys.setenv(LANG = "ja_JP.UTF-8")

する。

 

スピログラフ--地球と金星

Spirograph pattern: Venus x Earth
https://jp.mathworks.com/matlabcentral/communitycontests/contests/4/entries/3621

にある。

西洋占星術のページでは「地球と金星を結ぶ線分が薔薇を描く」と神秘的な取り扱いをしているところもある。

が，まあ，要するにスピログラフ Spirograph である。

Julia で書いて（描いて）みた。

using Plots
function roseofplanet(radius1, period1, radius2, period2; step=4, n=div(lcm(period1, period2), step), w=400, h=400)
    function circle(r, color=:green)
        θ = range(0, 360, length=1000)
        plot!(r .* cosd.(θ), r .* sind.(θ), color=color, lwd=0.4)
    end
    pyplot(grid=false, aspect_ratio=1, size=(w, h),
        xlabel="", ylabel="", label="", alpha=0.2)
    θ1 = range(0, step=360/period1*step, length=n)
    θ2 = range(0, step=360/period2*step, length=n)

    x1 = radius1 .* cosd.(θ1);
    y1 = radius1 .* sind.(θ1);
    x2 = radius2 .* cosd.(θ2);
    y2 = radius2 .* sind.(θ2);

    plot(showaxis=false);
    for (xi, yi, xj, yj) in zip(x1, y1, x2, y2)
        plot!([xi, xj], [yi, yj], color=:pink, lwd=0.2)
    end
    circle(radius1, :cadetblue1)
    circle(radius2, :yellow)
    plot!()
    savefig("fig.pdf")
end

引数

radius1, period1　片方の惑星の長半径と公転周期
radius2, period2　もう一方の惑星の長半径と公転周期
step=4　何日おきに描くか
n=div(lcm(period1, period2), step)　全描画日数
w=400, h=400　画像のピクセルサイズ

スピログラフの肝は，2個の惑星の公転周期が素であってはあまり良くないということだ。

例に挙げられている地球と金星の場合だと，実際の 365.2 日と 224.7 日を 364, 224 とするときれいな図になる。

rVenus, vVenus = 1.08, 224; # 1.082, 224.7
rEarth, vEarth = 1.50, 364; #1.495, 365.2
rMercury, vMercury = 0.579, 88; # 1.082, 224.7
rMars, vMars = 2.279, 690; # 2.279, 687

roseofplanet(rEarth, vEarth, rVenus, vVenus, step=4)

step=1 にすると，より本物らしく見える？
roseofplanet(rEarth, vEarth, rVenus, vVenus, step=1)

途中まで描画するとカージオイドが見える。
roseofplanet(rEarth, vEarth, rVenus, vVenus, step=1, n=600)

単なるスピログラフなので，金星と水星，火星と地球でも似たようなものができる。薔薇というより菊に近いか。
roseofplanet(rVenus, vVenus, rMercury, vMercury, step=2)

roseofplanet(rMars, vMars, rEarth, vEarth, step=50)

ネフロイドも見える。
roseofplanet(1, 100, 0.8, 300, step=1)

 

 

 

 

Google Colaboratory で Julia を使う（決定版！）

前報の方法はいささかダサいので，以下を決定版として推します。

Julia_Colab_Notebook_Template.ipynb - Colaboratory - Google
https://colab.research.google.com/github/ageron/julia_notebooks/blob/master/Julia_Colab_Notebook_Template.ipynb

上のリンクをクリックすると，Google Colaboratory で Julia_Colab_Notebook_Template.ipynb を開いた状態になる。

最初のテキストブロックに説明が書いてある。

要点は，

0. 上の方に

JULIA_VERSION="1.7.1" # any version ≥ 0.7.0

JULIA_PACKAGES="IJulia BenchmarkTools Plots"

があるが，必要なバージョンを指定することと，余計な時間がかかるので，BenchmarkTools, Plots は削除しておく

JULIA_VERSION="1.7.2" # any version ≥ 0.7.0

JULIA_PACKAGES="IJulia"

1. メニューバーの「ファイル」から「 ダウンロード」 --> 「.ipynb をダウンロード」により Julia_Colab_Notebook_Template.ipynb が保存される

2. メニューバーの「ファイル」から，「ノートブックをアップロード」でダウンロードした Julia_Colab_Notebook_Template.ipynb をアップロードする（バージョンが変更されていることを確認）

3. 次のコードブロックを実行する（ちょっと時間がかかる。50秒弱。）

4. ⌘+R で，Colaboratory のページをリロードする <-- これが最重要（これをしないと Julia カーネルにならない）

5. versioninfo() が書かれているコードブロックを実行する。

以下のようなものが出力されれば OK!! Julia が使える。

Julia Version 1.7.2
Commit bf53498635 (2022-02-06 15:21 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Xeon(R) CPU @ 2.20GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, broadwell)
Environment:
  JULIA_NUM_THREADS = 2
  

Google Colaboratory で Julia を使う（見なくていい）

決定版を見よ！
https://blog.goo.ne.jp/r-de-r/e/2c407c3a8cc68c4cd7986e4e5efa3069

 

 

 

 

この方法はちょっとダサい。

Google Colaboratory は Python しか使えないと思っている人が多いが，Julia も使える。

前報の R  を使うよりはちょっと手数が多い。

https://qiita.com/cometscome_phys/items/1ba6ec181bb0fe1b35d5

を要約すれば，

以下のような内容の雛形ファイル（拡張子 .ipynb）を作っておく（コピペして保存しておく）。

Julia のバージョンにより

      "name": "julia-1.7",
      "display_name": "julia-1.7"

      "!curl -sSL \"https://julialang-s3.julialang.org/bin/linux/x64/1.7/julia-1.7.2-linux-x86_64.tar.gz\" -o julia.tar.gz\n",

の該当 3 箇所は適切に変更する（例に上げているのは Julia 1.7.2 の場合）。

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "test.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "julia-1.7",
      "display_name": "julia-1.7"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "id": "CN08w8uiICa8",
        "colab_type": "code",
        "colab": {}
      },

      "source": [
        "!curl -sSL \"https://julialang-s3.julialang.org/bin/linux/x64/1.7/julia-1.7.2-linux-x86_64.tar.gz\" -o julia.tar.gz\n",
        "!tar -xzf julia.tar.gz -C /usr --strip-components 1\n",
        "!rm -rf julia.tar.gz*\n",
        "!julia -e 'using Pkg; pkg\"add IJulia\"'"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "sGOaN98-IZxj",
        "colab_type": "text"
      },
      "source": [
        "上記をshift+enterで実行後、「ランタイム」->「ランタイムのタイプを変更」->「保存」を押してください。"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ls5Cl0hsIndc",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "versioninfo()"
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}

このファイルを 「ファイル --> ノートブックのアップロード」 でアップロードし，書かれているように最初のコードブロックを実行する。

Julia をインストールスので，ちょっと時間がかかるが，辛抱だ。

インストールが終わったら，

「ランタイム --> ランタイムのタイプを変更 --> 保存」をすることを忘れないように。

その後，2 番目のコードブロック versioninfo() を実行してみよう。

Julia Version 1.7.2
Commit bf53498635 (2022-02-06 15:21 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Xeon(R) CPU @ 2.20GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, broadwell)

などと書かれれば，正常にインストールできていることがわかる。

以後，何でも書けばよい。

f(x) = 2x
f(8) |> println
f(1) |> println

16
2

 

Google Colaboratory で R を使う

Google Colaboratory は Python しか使えないと思っている人が多いが，R も使える。

以下のような内容の雛形ファイル（拡張子 .ipynb）を作り（コピペして保存しておく），

雛形ファイル.ipynb

{
  "cells": [],
  "metadata": {
    "colab": {
      "collapsed_sections": [],
      "name": "Title",
      "toc_visible": true,
      "provenance": []
    },
    "kernelspec": {
      "display_name": "R",
      "name": "ir"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}

「ファイル --> ノートブックをアップロード」で作っておいた雛形ファイルを指定して読み込む。

一度だけ，うまく行ったかの確認をするなら，

「ランタイム --> ランタイムのタイプを変更」で表示されるノートブックの設定ウインドウで，ランタイムのタイプが R になっていれば OK。

これで，R が使える。保存が必要なら，ファイル名を変えて「ファイル --> 名前の変更」から保存

「ファイル --> 保存」は Google Drive への保存になる。

「ファイル --> ダウンロード --> .ipynb をダウンロード」 で，ローカルドライブへの保存。

 

0/1 データの集計

各列が 0/1 の二値データ，グループ変数別，列ごとに 1 の個数を集計する。

データの生成と保存

サンプルサイズ n = 20000，列数 m = 20, グループ変数 二値データ

using RCall

R"""
set.seed(123)
n = 20000
m = 20
df = as.data.frame(matrix(sample(0:1, n*m, replace=TRUE), n, m),
    row.names=paste0("V", 1:m))
df$ID = 1:n
df$gender = factor(sample(c("male", "female"), n, replace=TRUE),
                    levels=c("male", "female"))
write.csv(df, "multianswer.csv", row.names=FALSE)
"""

tidyverse では， df[*, *] のような記述をひどく忌み嫌うとかで，

20  変数の集計は，20行書かなければならない（1000変数だったら，プログラムでプログラムを書く）。

R"""
library(tidyverse)
library(dplyr)
library(tidyr)

df = read_csv("multianswer.csv", show_col_types = FALSE)
system.time({
    cross <- df %>%
    group_by(gender) %>%
    summarise(
    V1 = sum(V1, na.rm = TRUE),
    V2 = sum(V2, na.rm = TRUE),
    V3 = sum(V3, na.rm = TRUE),
    V4 = sum(V4, na.rm = TRUE),
    V5 = sum(V5, na.rm = TRUE),
    V6 = sum(V6, na.rm = TRUE),
    V7 = sum(V7, na.rm = TRUE),
    V8 = sum(V8, na.rm = TRUE),
    V9 = sum(V9, na.rm = TRUE),
    V10 = sum(V10, na.rm = TRUE),
    V11 = sum(V11, na.rm = TRUE),
    V12 = sum(V12, na.rm = TRUE),
    V13 = sum(V13, na.rm = TRUE),
    V14 = sum(V14, na.rm = TRUE),
    V15 = sum(V15, na.rm = TRUE),
    V16 = sum(V16, na.rm = TRUE),
    V17 = sum(V17, na.rm = TRUE),
    V18 = sum(V18, na.rm = TRUE),
    V19 = sum(V19, na.rm = TRUE),
    V20 = sum(V20, na.rm = TRUE))
})
# cross
"""

RObject{VecSxp}

# A tibble: 2 x 21
  gender    V1    V2    V3    V4    V5    V6    V7    V8    V9   V10   V11   V12
1 female  4888  4954  4894  4940  4931  4964  4915  4974  4952  4978  5012  5017
2 male    5054  5042  5079  5065  4992  5028  5105  5021  5095  5011  5084  5109
# ... with 8 more variables: V13 , V14 , V15 , V16 ,
#   V17 , V18 , V19 , V20

RObject{RealSxp}
   user  system elapsed 
  0.004   0.001   0.005 

base-R で書けば，至極簡単。何万変数あったって，なんてことないね。

R"""
df = read.csv("multianswer.csv")
m = 20
system.time({i = m + 2
j = 1:m
ans <- sapply(j, function(k) table(df[, i], df[, k]))[3:4, ]
# ans

})
"""

RObject{IntSxp}
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 4888 4954 4894 4940 4931 4964 4915 4974 4952  4978  5012  5017  4901  5000
[2,] 5054 5042 5079 5065 4992 5028 5105 5021 5095  5011  5084  5109  5008  4927
     [,15] [,16] [,17] [,18] [,19] [,20]
[1,]  4856  4971  4936  4920  4987  4933
[2,]  5086  4991  5068  4986  4982  5095

RObject{RealSxp}
   user  system elapsed 
  0.027   0.001   0.028 

まあ，Julia だと 25 倍速い。とは言っても，1 秒未満なのだから競争しても意味ないが。

using CSV, DataFrames

function multianswer(df, m)
    gdf = groupby(df, :gender);
    vcat(sum(Matrix(gdf[2][!, 1:m]), dims=1), sum(Matrix(gdf[1][!, 1:m]), dims=1))
end

df = CSV.read("multianswer.csv", DataFrame);

@time multianswer(df, 20)

0.001102 seconds (413 allocations: 3.527 MiB)

2×20 Matrix{Int64}:
 4888  4954  4894  4940  4931  4964  …  4856  4971  4936  4920  4987  4933
 5054  5042  5079  5065  4992  5028     5086  4991  5068  4986  4982  5095

Julia のパイプ演算子 |>

Julia のパイプ演算子は引数が 1 個の関数を，一度に 1 個の関数のみをサポートする。

Pipe パッケージを使うことにより，base-R の |> と同じように，'_' をプレースホルダーとして使うことができるようになる。

前の記事 R: magrittr の %>% と base-R の |> で挙動が異なる
と同じ計算 (((2+3)-4)*5/6)^7 をやるためには，まず二項演算関数を定義する。

julia> add(x, y) = x + y 
julia> mul(x, y) = x * y
julia> sub(x, y) = x - y
julia> div(x, y) = x / y
julia> pow(x, y) = x ^ y

然る後，以下のようにする。

julia> using Pipe
julia> @pipe 2 |> add(_, 3) |> sub(_, 4) |> mul(_, 5) |> div(_, 6) |> pow(_, 7) |> println
0.2790816472336535

最後の println は括弧をつけない。付けるとすれば println(_) とする必要がある。

うん。あまりかっこよくないね。パイプを使うためにそれに使える関数を書くなどとは，本末転倒か?

R: magrittr の %>% と base-R の |> で挙動が異なる

以下の実行結果はなんの不審もないだろう。

> (((2+3)-4)*5/6)^7
[1] 0.2790816

magrittr パッケージの %>% を使って（無理やり計算してみる）

> library(magrittr)
> 2 %>% +(3)
[1] 5
> 2 %>% -(3)
[1] -1
> # 2 %>% *(3) # これはエラーになる。
> 2 %>% "*"(3)
[1] 6
> # 2 %>% /(3) # これはエラーになる。
> 2 %>% "/"(3)
[1] 0.6666667
# 2 %>% ^(3) # これはエラーになる。
> 2 %>% "^"(3)
[1] 8

まとめると
> 2 %>% +3%>% -4 %>% "*"(5) %>% "/"(6) %>% "^"(7)
[1] 0.2790816

base-R では，基本的に四則演算子はダブルクオートで区切る，
かつ第1引数は名前はなんでもよいが形式的に ' for=_' とでもする
> 2 |> "+"(a=_, 3)
[1] 5
> 2 |> "-"(any=_, 3)
[1] -1
> 2 |> "*"(foo=_, 3)
[1] 6
> 2 |> "/"(longlonglonglonganything=_, 3)
[1] 0.6666667
> 2 |> "^"(xxxxx=_, 3)
[1] 8

まとめると
> 2 |> "+"(d=_, 3) |> "-"(d=_, 4) |> "*"(d=_, 5) |> "/"(d=_, 6) |> "^"(d=_, 7)
[1] 0.2790816

同じ結果を得るにしても，かなり指定法が異なる。

なお，Julia の |> については，調査結果が得られていない。どうしたもんか...