行列とベクトルの積

twitter はやらないのでこちらでコメント

R の行列は列優先なので，

x = matrix(1, 3, 3)
y = 1:3

のときに，x*y でよい場合と，t(t(x)*y) が欲しい場合がある。

> x*y
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    2    2    2
[3,]    3    3    3

> t(t(x)*y)
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    1    2    3
[3,]    1    2    3

後者の場合に t(t(*)) がダサいと思うので，何かよい方法がないのか...ということだが。

x* matrix(rep(y, 3), 3, byrow=TRUE)

のようにするのも一法だ。

> x* matrix(rep(y, 3), 3, byrow=TRUE)
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    1    2    3
[3,]    1    2    3

しかし，y を行方向に水増しした行列をつくって掛け算するということなので，無駄であることは明らか。x の行数が小さい場合にはやや有利なのかもしれない。

> x = matrix(1, 3, 3)
> y = 1:3
> system.time(replicate(1e4, t(t(x)*y)))
   ユーザ   システム       経過  
     0.215      0.008      0.242
> system.time(replicate(1e4, x* matrix(rep(y, 3), 3, byrow=TRUE)))
   ユーザ   システム       経過  
     0.050      0.000      0.049

> x = matrix(1, 300, 3)
> y = 1:3
> system.time(replicate(1e4, t(t(x)*y)))
   ユーザ   システム       経過  
     0.360      0.040      0.406
> system.time(replicate(1e4, x* matrix(rep(y, 300), 300, byrow=TRUE)))
   ユーザ   システム       経過  
     0.321      0.027      0.348

> x = matrix(1, 3000, 3)
> y = 1:3
> system.time(replicate(1e4, t(t(x)*y)))
   ユーザ   システム       経過  
     2.315      0.774      3.074
> system.time(replicate(1e4, x* matrix(rep(y, 3000), 3000, byrow=TRUE)))
   ユーザ   システム       経過  
     2.951      0.909      3.844

10 を作る（その 3）

3478 で 10 を作れと...

難易度，星二つでも，最難関でも，プログラム書いておいたら，そんなの関係ね～～

(3 - (7 / 4)) * 8
8 * (3 - (7 / 4))

でもなあ。今回は，今までにない条件が付いている。

ルール１：数字を並び替えることは出来ない。

これでは，解はないのではないか？

2015/01/13 19:00 時点で，以下の様な修正が記載されていた。いい加減にしろ。

 

（注：ルール１を誤って並び替えることは出来ないとしていました。申し訳ありません。）

もっとも，

マシュー：「ぜひ挑戦してみたいな。これができたらいっしょにデートに行かないか？」
リンカ　：「そうね。天才とならデートに行きたいわ。」

ということなので，例によって，はなから相手にされていないのかもね。

ちなみに，「●●を並び替える」は誤った日本語であり，「●●を並べ替える」でなくてはならない。

「●●を並び替える」，「●●の並び替え」と誤記しているページはいっぱいある。なげかわしい。

そもそも，『「並び替え」と「並べ替え」の違いは何か』などという質問まであるのが，なげかわしい。

規則性のある文字列出力（3）

R は，短くする方策がまだ見つからないので，現実逃避で AWK で書いてみる

BEGIN{ } は除くようなので，41文字解ができた（同率首位のようだ）

暗号化して掲示しておく

Grsvo(s<cai)zBsxDp"%m",s%cg+(s++%i?jh:gf)

締め切りが過ぎたので，平文化する

while(i<208)printf\"%c\",i%26+(i++%8?97:65)

規則性のある文字列出力（2）

37 文字解があるようだ。驚異的だなあ。

45 文字解（暗号化した）

mkD(BkG3yMrkB(kC.BkG(gf:ja+dc*(a:cah%%i>a))))

締め切りが過ぎたので，平文化する

cat(rawToChar(as.raw(65:90+32*(0:207%%8>0))))

規則性のある文字列出力

AbcdefghIjklmnopQrstuvwxYzabcdefGhijklmnOpqrstuvWxyzabcdEfghijklMnopqrstUvwxyzabCdefghijKlmnopqrStuvwxyzAbcdefghIjklmnopQrstuvwxYzabcdefGhijklmnOpqrstuvWxyzabcdEfghijklMnopqrstUvwxyzabCdefghijKlmnopqrStuvwxyz

を出力する

暗号化した解

mkD(zkCDo(spovCo(a:cah%%i,voDDoBC,VO33O12),myvvkzCo=""))

回答締め切りが過ぎたので，平文化

cat(paste(ifelse(0:207%%8,letters,LETTERS),collapse=""))

strptime が面白い

strptime が面白い

> strptime("2014年12月27日 14時27分34秒", "%Y年%m月%d日 %H時%M分%S秒")
[1] "2014-12-27 14:27:34 JST"

> strptime("2014/12/27 14-27-34", "%Y/%m/%d %H-%M-%S")
[1] "2014-12-27 14:27:34 JST"

> strptime("2014.12.27 14.27.34", "%Y.%m.%d %H.%M.%S")
[1] "2014-12-27 14:27:34 JST"

> strptime("14-27-34 2014/12/27", "%H-%M-%S %Y/%m/%d")
[1] "2014-12-27 14:27:34 JST"

> strptime("14-27-34 2014/12/27", "%H-%M-%S %Y/%m/%d")
[1] "2014-12-27 14:27:34 JST"

> strptime("14-27 2014/12/27", "%H-%M %Y/%m/%d")
[1] "2014-12-27 14:27:00 JST"

> strptime("2014/12/27", "%Y/%m/%d ")
[1] "2014-12-27 JST"

> strptime("14-27-34", "%H-%M-%S")
[1] "2015-01-08 14:27:34 JST"

適合度の検定について（2）

「適合度の検定について」に対し，中澤さんからコメントをいただいた。

> vcd パッケージの goodfit が便利です

はじめて使ってみました。

> goodfit(as.table(tab1))

Observed and fitted values for poisson distribution
with parameters estimated by `ML'

 count observed       fitted
     0        0 2.150092e-03
     1        0 2.805870e-02
     2        0 1.830830e-01
     3        1 7.964111e-01
     4        1 2.598291e+00
     5        3 6.781540e+00
     6       12 1.474985e+01
     7       27 2.749793e+01
     8       42 4.485600e+01
     9       63 6.504121e+01
    10       80 8.487877e+01
    11      108 1.006971e+02
    12      129 1.095081e+02
    13      108 1.099293e+02
    14      106 1.024698e+02
    15       87 8.914871e+01
    16       79 7.271192e+01
    17       45 5.581709e+01
    18       42 4.046739e+01
    19       31 2.779471e+01
    20       18 1.813605e+01
    21        8 1.127026e+01
    22        5 6.685312e+00
    23        1 3.793188e+00
    24        1 2.062546e+00
    25        2 1.076649e+00
    26        1 5.403950e-01

> summary(goodfit(as.table(tab1)))

     Goodness-of-fit test for poisson distribution

                      X^2 df  P(> X^2)
Likelihood Ratio 19.73707 22 0.5994826

fitted の合計が 999.5218 となるのがちょっといやな感じです（x >= 27 に対する確率が考慮されていない）。

goodfit がやっていることは，以下の通り（冗長な部分を削除）。

x = 3:26
tab1 = c(1, 1, 3, 12, 27, 42, 63, 80, 108, 129, 108, 106, 87, 79, 45, 42, 31, 18, 8, 5, 1, 1, 2, 1)
names(tab1) = x
lambda = sum(x * tab1) / sum(tab1)
p = dpois(x, lambda)
cbind(x, tab1, p, p*sum(tab1))
ln = function(o, e) sum(ifelse(o == 0, 0, o*log(o/e)))
( G2 = sum(ln(tab1, p*sum(tab1)))*2 )
( df = length(x) - 1 -1 )
( p.value = pchisq(G2, df, lower.tail=FALSE) )

結果

> ( G2 = sum(ln(tab1, p*sum(tab1)))*2 )
[1] 19.73707
> ( df = length(x) - 1 -1 )
[1] 22
> ( p.value = pchisq(G2, df, lower.tail=FALSE)　)
[1] 0.5994826

G2（G スクエアー）を使う検定（独立性の検定でも，G2 を使う検定が定義されている）。

この検定では，観察値が 0 のセルは検定統計量に反映されないので，観察値が 0 のセルがあっても良いし，「期待値が小さいセルについての考慮」も不要ということにはなっている。

適合度の検定について

「交代再生過程のある時刻における状態がポアソン分布に従うか否かシミュレーションしてみた.(非公開にされました)」において，

> 検定の仕方まちがえてるのだろうか？

ということですが。そのとおり，残念ながら，間違えています。

著者が最初にやった chisq.test(tab1) は，一様分布に従うかどうかの検定です。一様とはとてもいえないので，帰無仮説は棄却されてしまいますね。

やるべき検定は，適合度検定です。chisq.test 関数を使います。

一様分布に従うかどうかも，適合度検定ですが，chisq.test の引数に p というのがあり，分布が理論比 p に従うかどうかの適合度の検定を行うのです。n 個のカテゴリーがある場合の一様分布の適合度検定は p = 1/n がデフォルトになっているのです。

ですから，ポアソン分布などへの適合度検定の場合には p を明示的に指定してやらなければなりません。これが，一筋縄ではいかないことがあります。

まず，最初にポアソン分布に従うという帰無仮説のもとで適合度検定を行うわけですが，ポアソン定数が分かっているならそれを使えばよいのですが，もし標本から推定したポアソン定数（平均値）を使うならば，後述するカイ二乗分布の自由度が余分に 1 減ります。

確率変数（例の場合は 3 ～ 26）をとる確率は dpois で求まります（x1 = 3:26); dpois(x1, l1)）が，ポアソン分布の確率変数は 0 から ∞ ですが，0 近辺や 大きな値をとる確率は小さくなるので，小さい方は 0 から ある値までの確率を合計します（プールします）。同様にある値以上の確率も合計します。「ある値」は，期待値が 1 以上になるような値です。幸いなことに例題の場合は，最小と最大（3 と 26）になります。期待値が 1 以上になるようにするためにはもっと大きい，あるいは，小さい値までを含めるなければならないこともあるでしょう。この場合には，度数分布表もまとめる必要があるでしょう（そして，そのような場合には結果として自由度が余分に減少します）。
ちょっとややこしいですが，丁寧に説明してあるページなどを参照しましょう。

で，最終的には以下に示すデータフレームのような数値を使って検定を行います。

         x1 tab1       Dpois Expectation
0 thru 3  3    1 0.001009703    1.009703
4         4    1 0.002598291    2.598291
5         5    3 0.006781540    6.781540
6         6   12 0.014749849   14.749849
7         7   27 0.027497933   27.497933
8         8   42 0.044856004   44.856004
9         9   63 0.065041206   65.041206
10       10   80 0.084878773   84.878773
11       11  108 0.100697090  100.697090
12       12  129 0.109508086  109.508086
13       13  108 0.109929271  109.929271
14       14  106 0.102469784  102.469784
15       15   87 0.089148712   89.148712
16       16   79 0.072711919   72.711919
17       17   45 0.055817090   55.817090
18       18   42 0.040467391   40.467391
19       19   31 0.027794708   27.794708
20       20   18 0.018136047   18.136047
21       21    8 0.011270258   11.270258
22       22    5 0.006685312    6.685312
23       23    1 0.003793188    3.793188
24       24    1 0.002062546    2.062546
25       25    2 0.001076649    1.076649
26 over  26    1 0.001018650    1.018650

自由度は，カテゴリー数（階級数）マイナス 2 ですが，マイナス 1 は度数の合計がサンプルサイズになるという制約，もう一つのマイナス 1 は，標本からポアソン定数を決めたための制約です。

ということで，結果は，カイ二乗値=16.18529，自由度=22，P値=0.8065844ということになり，「ポアソン分布に従っていないとはいえない」という結論になります。

なお，二項分布への適合度検定も同じように行えますが，その場合も，帰無仮説は棄却できないようです。

x1 = 3:26
tab1 = c(1, 1, 3, 12, 27, 42, 63, 80, 108, 129, 108, 106, 87, 79, 45, 42, 31, 18, 8, 5, 1, 1, 2, 1)
n = sum(tab1)

# chisq.test(d$freq) # これは一様分布の検定なので間違い

( l1 = sum(x1 * tab1) / n ) # lambda
Dpois = dpois(x1, l1)
Dpois[1] = sum(dpois(0:3, l1)) # = ppois(3, l1)
Dpois[24] = 1-sum(Dpois[1:23]) # = ppois(25, l1, lower.tail=FALSE)
Expectation = Dpois * n
d = data.frame(x1, tab1, Dpois, Expectation)
rownames(d) = c("0 thru 3", 4:25, "26 over")
d
( chisq = sum((tab1-Expectation)^2/Expectation) )
( df = length(tab1) - 1 - 1 )
( p.value = pchisq(chisq, df, lower.tail=FALSE) )

結果

> ( chisq = sum((tab1-Expectation)^2/Expectation) )
[1] 16.18529
> ( df = length(tab1) - 1 - 1 )
[1] 22
> ( p.value = pchisq(chisq, df, lower.tail=FALSE) )
[1] 0.8065844

R で移動平均

「R で移動平均を求める--R プログラミングの小ネタ」で，

> Rには移動平均そのものずばりを求める関数はないようです。

とあるが，latticeExtra の simpleSmoothTs というのがある。width はデフォルトでは NROW(x) %/%10 + 1 なので，自分で決めた方がよいようだが。

x = c(67, 59.2, 48.4, 56.5, 66.9, 79.1, 91.8, 87.6, 52.2, 66.4, 99.3,
159.3, 74.3, 57.3, 38.5, 42.7, 46.2, 42, 24.7, 26.4, 43.4, 38.5,
44, 31.8, 12.5, 41.9, 31.3, 44.2, 26.6, 23.2, 18, 20.8, 33, 9.4,
16.6, 11.1, 15.2, 10.9, 9.6, 8.3, 14.4, 18.3, 8, 8.6, 8.1, 6.9,
10, 10.4, 10.8, 19.1, 7.2, 2.2, 1, 4.2, 3.2, 0, 5.6, 3.9, 0.9,
0, 3, 7.4, 3.5, 12.8, 20.8, 5.9, 3, 7.2, 6.8, 6.6, 7.9, 13.3,
27.9, 26.5, 35.1, 19, 36.9, 66, 43.1, 42.7, 90.9, 108.8, 111.4,
94, 75.6, 75.1, 77.9, 78.4, 59.3, 68, 46.4, 76.2, 71.4, 72.2,
79.2, 112.3, 98.1, 117.5, 70.3, 94.2, 99.7, 116.7, 134.4, 136.7,
100.8, 117, 80.7, 88.5, 97.7, 123.6, 97.8, 119.3, 143.9, 111,
175.7, 143.1, 128.3, 130, 140.7, 105.7, 90.1, 83.6, 67.9, 123.5,
124.5, 112.8, 110, 121.5, 135.4, 106.1, 77.8, 57.2, 54.1, 63.1,
84, 70.7, 68, 89.1, 73.4, 110.2, 74.1, 61.1, 67.3, 70.7, 46.2,
43.4, 46.8, 31.5, 31.4, 66.3, 79.4, 70.4, 77.4, 50.9, 43.8, 39.3,
38, 48.7, 56.2, 68.8, 26.1, 11.2, 13.3, 21.8, 17, 19.2, 31.2,
24.8, 13.1, 40.9, 27.6, 32.3, 11.7, 9.1, 14.8, 16.2, 11.7, 6,
8.2, 14.4, 1.4, 0.9, 9.8, 0, 1.7, 4.1, 2, 14, 10.3, 12.8, 9.7,
23.7, 17.1, 6.9, 19.4, 33.3, 29.9, 35.1, 30.8, 25.2)

library(latticeExtra)
plot(x, type="l", lty=3)
lines(simpleSmoothTs(x), col=2, lwd=2)
lines(simpleSmoothTs(x, width=11), col=4, lwd=2)

その他にも，いくつかあるが...それぞれの関数の仕様の違いに注意すべし。

m を整数としたとき，パッケージ TTR の SMA(x, m+1) は simpleSmoothTs(x, width=m, sides=1) と同じ結果になる。

標準パッケージ stats の filter(x, rep(1, m))/m は simpleSmoothTs(x, width=m) と同じ結果になる。

2つのベクトル相互で共通しない要素のリストアップ

a = c("foo", "bar", "baz", "abc", "de")
b = c("hoge", "hogehoge", "abc", "de")

%in% と "[" を使うと

> a[! a %in% b]
[1] "foo" "bar" "baz"

> b[! b %in% a]
[1] "hoge"     "hogehoge"

しかし，setdiff を使うのがモア・ベター？？

> setdiff(a, b)
[1] "foo" "bar" "baz"

> setdiff(b, a)
[1] "hoge"     "hogehoge"

円周率の連分数展開

「円周率の連分数展開」であるが，22 段階まで計算すれば有効桁数 16 の近似値が得られる。

n = 22
y = 2*1:n-1
x = c(4, (1:(n-1))^2)
Pi = 0
for (i in n:1) {
    Pi = x[i]/(y[i]+Pi)
}

> Pi
[1] 3.141592653589793
> options(digits=16)
> pi
[1] 3.141592653589793

属性相関係数

「効果量としてのファイ係数とクラメールのV」だが

引用している「研究論文における効果量の報告のために」において

> 2×2 の分割表（クロス表）の場合には，相関係数の一種であるφ（ファイ）係数を用い，それ以外の場合は，Cramer's V を利用する。

不完全な記述だなあ。

2×2分割表の場合には「相関係数の一種」ではなく，ピアソンの積率相関係数そのもの（ただし，定義の上ではピアソンの積率相関係数の絶対値。方向性を自分で付ければ，相関係数そのもの）

φ係数は2×2分割表以外の場合にも定義できる。ただ，2×2分割表以外の場合には1以上の値を取るので，相関係数の場合と同じように0～1の範囲になるように制限した（正規化した）ものがクラメールのV