上級者のための R および RjpWiki に関する質問コーナー
注意:新規記事用の入力欄は以下の目次の直後にあります。各記事への追加コメントの入力欄は、各記事の後にあります。質問が長くなる場合は、まず出だしだけ投稿し、次に編集ボタンでそれを修正、追加するのが便利です。
ここでは、甘えは一切許されません。必ずサマリーページ(もしくは本家への投稿)を作る必要があります。 初心者の方でも投稿は構いませんが
恐がる必要はありません。答える以上、それがまとまった情報になることを望んでいるだけです。
青木繁伸 (2008-12-07 (日) 12:43:40)
cor 関数については,遙か昔に「欠損値を含むデータにおいて誤った結果を返す」というバグがあった。R 2.8.0 から use 引数で指定できるものが従来の3つから以下の5つに増えた。
"everything", "all.obs", "complete.obs", "na.or.complete", or "pairwise.complete.obs".
この追加にさいして,またバグが持ち込まれた。
問題が生じるのは,確認した範囲では use="complete.obs" で method="spearman" のときである。
現状は以下の通り。wrong result ----- BUG in R 2.8.0 > cor(airquality[1:4], use="complete.obs", method="spearman") Ozone Solar.R Wind Temp Ozone 1.0000000 0.35846998 -0.60797155 0.7807127 Solar.R 0.3584700 1.00000000 -0.07022037 0.2157419 Wind -0.6079715 -0.07022037 1.00000000 -0.5034014 Temp 0.7807127 0.21574193 -0.50340142 1.0000000この結果は,間違いである。use="complete.obs" は要するに,NA を1つでも含む行は削除した上で完全なデータについて相関係数を求めるということであるから,以下のようにして得られるのと同じ結果でなければならない。
> cor(na.omit(airquality[1:4]), use="complete.obs", method="spearman") Ozone Solar.R Wind Temp Ozone 1.0000000 0.34818647 -0.60513642 0.7729319 Solar.R 0.3481865 1.00000000 -0.06169636 0.2095369 Wind -0.6051364 -0.06169636 1.00000000 -0.4993228 Temp 0.7729319 0.20953692 -0.49932278 1.0000000食い違いのある結果のいずれが正しいかを検証するには,たとえば,Ozone と Temp のスピアマンの順位相関係数をとるために,以下のようにした結果と比べればよい。sample estimates は 0.772932 である。
> a <- subset(airquality[1:4], complete.cases(airquality[1:4])) > cor.test(a[,1], a[,4], method="spearman") Spearman's rank correlation rho data: a[, 1] and a[, 4] S = 51753.35, p-value < 2.2e-16 alternative hypothesis: true rho is not equal to 0 sample estimates: rho 0.772932さて,では現在の cor 関数のどこがまずく,どのように直せばよいのかについて,検討したところ,以下の部分をなおせば正しい結果が得られるようになることを確認した。修正箇所を示すために,少し広めの範囲を引用しておく。
else if (na.method != 3L) { Rank <- function(u) { if (length(u) == 0L) u else if (is.matrix(u)) apply(u, 2L, rank, na.last = "keep") else rank(u, na.last = "keep") } ### begin ### # for 'use = "complete.obs"', add two lines below if (na.method == 2L) x <- na.omit(x) ### end ### x <- Rank(x) if (!is.null(y)) y <- Rank(y) .Internal(cor(x, y, na.method, method == "kendall")) }この修正後は,以下のように正しい結果が得られる。
> cor(airquality[1:4], use="complete.obs", method="spearman") Ozone Solar.R Wind Temp Ozone 1.0000000 0.34818647 -0.60513642 0.7729319 Solar.R 0.3481865 1.00000000 -0.06169636 0.2095369 Wind -0.6051364 -0.06169636 1.00000000 -0.4993228 Temp 0.7729319 0.20953692 -0.49932278 1.0000000なお,method="kendall" のときには問題が生じないようである(詳細には検討していない) また,use="complete.obs" 以外についても,未検討である。
なるほど (2006-12-11 (月) 10:14:44)
初めて投稿させていただきます。よろしくお願いします。
RのパッケージであるSSPIRを用いて, カルマンフィルタと拡張カルマンフィルタを行いたいのです。demoには4つの事例があります。カルマンフィルタと拡張カルマンフィルタに相当するものがどれなのかよくわかりません。
ご教示下さい。
hina (2005-06-03 (金) 21:55:59)
1変量の時系列データにカルマンフィルタをかけるために、KalmanLikeという関数を使いたいのですが、関数の引数のmodに与えられるlistのそれぞれのデータの型が分かりません。
helpや論文を調べてみたのですが、具体的な例が見つかりません。
ご教示下さい。
'a <- T a + R e', e ~ N(0, kappa Q) システムモデル 'y = Z'a + eta', eta ~ N(0, kappa h) 観測モデル従ってmodではこれらのモデルの係数類(h以外はベクトルか行列)を指定することになります。
'T' 状態遷移行列(p×p) 'Z' 観測係数ベクトル(p×1) 'h' 観測ノイズの分散 'V' RQR'(p×p) 'a' 状態の初期値(p×1) a[0|0] 'P' 状態の初期分散共分散行列(p×p) P[0|0] 'Pn' P[t|t-1] フィルタリングで使用されることになるPです。nitで指定された時点までこの値が固定的に使われます。 nit=0の場合は最初の1回だけ使われて、以降は普通に更新されていきます。余談ですがこのモデルではZが固定(時変でない)なのがかなり適用範囲を狭めていると思います。これだと単純な時変係数の回帰分析もできません。それで私はZを行列で受け付けるように改造したものを使っています。 -- makoto 2006-01-16 (月) 01:18:26
Uwamy@mpc (2005-06-02 (木) 11:33:35)
先日JIS,ISO式四捨五入を読み参考にしました。
しかしながら 12.95を小数第一位で丸めたいと思っていますが、
round(12.95,1)の結果は12.9となります。(正解は13.0のはず)
これってRのBUG?
なおR2.1.0を使いました
round(12.95+1e-15, 1)を試してみましょう。 -- 青木繁伸 2005-06-02 (木) 12:15:46
sum(rep(12.95,100)) options(digits=22) sum(rep(12.95,100))などを実行してみましょう.近似値自体の丸め処理もIEEEに基づいてますので roundされたものにroundして値がおかしいと言われても困る訳でして.
round(12.95*10^2,1-2)みたいな方法でもいいのかも. あと、ここはパッチを作成するなりして本家に報告せねばならぬと言う 恐ろしい場所なのですが...:-P -- なかま 2005-06-02 (木) 12:52:35
R初心者 (2004-09-27 (月) 22:40:36)
ExcelのデータをCSVに変換して読み込みを行った所,"hoge","foo","bar\nboo" # \nは実際の改行コード 1,2,3 4,5,6
> read.csv("a.csv") hoge foo bar.boo 1 boo\n1,2,3\n4,5,6\n NA NA
などとなってしまいます。セル内改行を除けば正常に読み込めます。
read.csv->read.table->scanと呼ばれていますが,Quoteに対する処理が見当たりませんでした。
Excelにおけるセル内改行してあるようなデータを正確に読む方法はないのでしょうか。