上級者のための R および RjpWiki に関する質問コーナー

注意:新規記事用の入力欄は以下の目次の直後にあります。各記事への追加コメントの入力欄は、各記事の後にあります。質問が長くなる場合は、まず出だしだけ投稿し、次に編集ボタンでそれを修正、追加するのが便利です。

ここでは、甘えは一切許されません。必ずサマリーページ(もしくは本家への投稿)を作る必要があります。 初心者の方でも投稿は構いませんが

恐がる必要はありません。答える以上、それがまとまった情報になることを望んでいるだけです。







R 2.8.0 の cor 関数での特定条件下のバグ

青木繁伸 (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" 以外についても,未検討である。

  • ずっとフォローしてなかったけど,R 2.13.0 では直っていた。 -- 青木繁伸 2011-05-31 (火) 14:44:54

SSPIR カルマンフィルタについて

なるほど? (2006-12-11 (月) 10:14:44)

初めて投稿させていただきます。よろしくお願いします。
RのパッケージであるSSPIRを用いて, カルマンフィルタと拡張カルマンフィルタを行いたいのです。demoには4つの事例があります。カルマンフィルタと拡張カルマンフィルタに相当するものがどれなのかよくわかりません。
ご教示下さい。

  • 改めて,上の方に書いてあることを読んでみましたが,ここって,すんごく,こわいところなんですね。認識してましたか? -- 2006-12-11 (月) 21:30:22

KalmanLike?の引数のmodについて

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,Z,h,R,Q,a)
    状態'a'の次元をpとすると、
    '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
  • help にはこの関数は StructTS 等の他の関数から呼ばれることを期待していると書かれています(だから example もない)が、そちらを参照されましたか。StructTS 自身(もしくはその example )コードの出力を点検されたらもう少し情報が得られるかも知れません。既に吟味済みでしたらごめんなさい。 -- つかったことない人? 2006-12-11 (月) 19:56:00

round関数 合ってますか?

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を使いました

  • バグなどではありません。
    「コンピュータ上に保持される実数は(整数や小数部が2^(-n)の和である場合を除き,要するに2進数で正確に表される場合以外),近似値である」
    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

Excelのデータ読み込みでエラー(例です)

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におけるセル内改行してあるようなデータを正確に読む方法はないのでしょうか。

  • scan.cにおけるデリミタの処理を適切に処理するようハックしてください。 -- R上級者? 2004-09-27 (月) 22:48:35
  • 週末ハックします。 -- R初心者? 2004-09-27 (月) 22:49:55



トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Google
WWW を検索 OKADAJP.ORG を検索
Last-modified: 2015-03-01 (日) 01:15:59 (1718d)