「統計学関連なんでもあり」の過去ログ--- 043

No.11314　いい関数を教えて下さい。　　【のざわ】　2009/11/23(Mon) 18:51

いつもお世話になっております。

よろしくお願いします。

Rで以下のようなSingular Spectrum Analysisのプログラムを作ったのですが，以下のような例では，きちんと一期先を返してきます。

> ssa(1:100)
[1] 101
> ssa(sin(1:100))
[1] 0.4520258
> sin(101)
[1] 0.4520258

しかし以下のようにすると，あさってなくらい吹き飛ばされてしまいます。

> ssa(Nile)
[1] 544.1812

たぶん，Nileがあちこちで不連続だから，吹き飛ばされるのかな？と思い以下のようにすると，今度は行列が計算されなくなってしまいます。

> ssa(spline(Nile))
Error in matrix(0, row, colmn) : invalid 'ncol' value (< 0)

何かこういう不連続なデータを旨く前処理して微分可能にして，行列も計算してくれるような尤度の高い補間関数があれば，お教え願えれば幸いです。

ざっくりRの中を見るとsplineくらいしか見当たりませんでした。

よろしくお願いします。

------------------------------------------------------------------

ssa <- function(data, row=10)
{
  d.len <- length(data)
  colmn <- d.len - row + 1
  size <- row*colmn
  X <- matrix(0, row, colmn)
  d.point <- 1
  for (L in 1:size) {
    X[L] <- data[d.point]
    if (L%%row == 0) {
      d.point <- d.point - row + 2
    } else {
      d.point <- d.point + 1
    }
  }
  C <- X%*%t(X)/colmn
  eigV <- eigen(C)$vectors
  V <- eigV[, 1:(row - 1)]
  L <- matrix(0, row, 1)
  L[row, 1] <- 1
  Q <- matrix(0, row, 1)
  Q[1:(row-1), 1] <- X[2:row, colmn]
  Xpre <- (( t(L)%*%V%*%t(V)%*%Q )/( 1- t(L)%*%V%*%t(V)%*%L))[1, 1]
  return(Xpre)
}

No.11315　Re: いい関数を教えて下さい。　　【のざわ】　2009/11/23(Mon) 19:18

少し補足説明します。

元々これは不連続があるものは，駄目なようです。
> ssa(tan(1:100))
[1] 74.55493
> tan(101)
[1] 0.5067526
tan()のように明らかに，切れ目がある奴は覚悟していましたが，Nileみたいなデータを旨く加工して，それとなく値を返して来て欲しいんですが。。

No.11316　Re: いい関数を教えて下さい。　　【青木繁伸】　2009/11/23(Mon) 20:02

解答ではないですが，ssa(spline(Nile)) がエラーになる理由は，spline が何を返すかをチェックするとわかるでしょう。ssa(spline(Nile)$y) とすれば，734.9737 が帰ってきますが，それがあなたのお気に召すかどうかはわかりません。

No.11317　Re: いい関数を教えて下さい。　　【のざわ】　2009/11/23(Mon) 20:13

いつも，どことなく間抜けな質問ばかりで済みません。spline関数の使い方を誤っていただけのようですね。。取り合えず暫く使ってみます。有難うございました。

No.11318　Re: いい関数を教えて下さい。　　【のざわ】　2009/11/23(Mon) 21:27

すみません。もうひとつ質問させてください。
こういうプログラムにデータを渡す場合，不連続を連続と見せたい場合の前処理は，他には，どんなものがあるでしょうか？
確かにsplineだと，おとなしくはなったのですが，高次な埋め込み次元ほど精度良く同定されたような結果になり，どうもイメージとは合わなくなってしまいました。

No.11319　Re: いい関数を教えて下さい。　　【のざわ】　2009/11/23(Mon) 22:57

どうも，やはりイメージと合わないので，こうやってみました。

> plot(Nile[1:100])
> lines(spline(Nile)$y)

ぜんぜん，splineは，fitしてる感じでは無いので，駄目なのか。と言う結論になりました。

No.11320　Re: いい関数を教えて下さい。　　【青木繁伸】　2009/11/23(Mon) 23:14

> ぜんぜん，splineは，fitしてる感じでは無いので，

あの～，length(spline(Nile)$y) をやってみて下さい。spline(Nile)$x がどのようになっているかも確かめて下さい。正しいグラフを描くには，以下のようにしないとダメですよ。（図はクリックすると原寸大に表示）
plot(Nile, type="l")
lines(spline(Nile), col=2)

No.11322　Re: いい関数を教えて下さい。　　【のざわ】　2009/11/24(Tue) 14:40

＞正しいグラフを描くには，以下のようにしないとダメですよ。

了解しました。
splineってxとyの両方が，各々幾らなのかを出力してきてるわけですね。暫くはsplineで色々やってみます。また何か結果が得られたら報告致します。
いつも，どうも有難うございました。