2006年08月18日

spec.ar

Rのspec.ar関数。なにかで修正しなきゃいけんかったんじゃけど、なぜ修正が必要だったか忘れた。

spec.ar <- function (x, n.freq, order = NULL, plot = TRUE, na.action = na.fail,
method = "yule-walker", ...)
{
if (!is.list(x)) {
series <- deparse(substitute(x))
x <- na.action(as.ts(x))
xfreq <- frequency(x)
nser <- NCOL(x)
x <- ar(x, is.null(order), order, na.action = na.action,
method = method)
}
else {
cn <- match(c("ar", "var.pred", "order"), names(x))
if (any(is.na(cn)))
stop("'x' must be a time series or an ar() fit")
series <- x\$series
xfreq <- x\$frequency
if (is.array(x\$ar))
nser <- dim(x\$ar)[2]
else nser <- 1
}
order <- x\$order
if (missing(n.freq))
n.freq <- 500
freq <- seq(0, 0.5, length = n.freq)
if (nser == 1) {
coh <- phase <- NULL
if (order >= 1) {
cs <- outer(freq, 1:order, function(x, y) cos(2 *
pi * x * y)) %*% x\$ar
sn <- outer(freq, 1:order, function(x, y) sin(2 *
pi * x * y)) %*% x\$ar
spec <- x\$var.pred/(xfreq * ((1 - cs)^2 + sn^2))
}
else spec <- rep(x\$var.pred/(xfreq), length(freq))
}
else .NotYetImplemented()
freq<-freq*xfreq
spg.out <- list(freq = freq, spec = spec, coh = coh, phase = phase,
n.used = nrow(x), series = series, method = paste("AR (",
order, ") spectrum ", sep = ""))
class(spg.out) <- "spec"
if (plot) {
plot(spg.out, ci = 0, ...)
return(invisible(spg.out))
}
else return(spg.out)
}



$の前に¥が見えたら消してくださいね(^_^;)
posted by metola at 14:35| Comment(0) | TrackBack(0) | 統計解析システムR | このブログの読者になる | 更新情報をチェックする
この記事へのコメント
コメントを書く
お名前: [必須入力]

メールアドレス:

ホームページアドレス:

コメント: [必須入力]

認証コード: [必須入力]


※画像の中の文字を半角で入力してください。

この記事へのトラックバック

Physical Review

Recent stories in Physical Review Focus

×

この広告は1年以上新しい記事の投稿がないブログに表示されております。