Como simular dados funcionais?

12

Estou tentando testar várias abordagens funcionais de análise de dados. Idealmente, eu gostaria de testar o painel de abordagens que tenho em dados funcionais simulados. Eu tentei gerar FD simulado usando uma abordagem baseada em um somar ruídos gaussianos (código abaixo), mas as curvas resultantes parecem muito resistentes em comparação com a coisa real .

Fiquei me perguntando se alguém tinha um ponteiro para funções / idéias para gerar dados funcionais simulados com aparência mais realista. Em particular, estes devem ser suaves. Sou completamente novo neste campo, portanto qualquer conselho é bem-vindo.

library("MASS")
library("caTools")
VCM<-function(cont,theta=0.99){
    Sigma<-matrix(rep(0,length(cont)^2),nrow=length(cont))
    for(i in 1:nrow(Sigma)){
        for (j in 1:ncol(Sigma)) Sigma[i,j]<-theta^(abs(cont[i]-cont[j]))
    }
    return(Sigma)
}


t1<-1:120
CVC<-runmean(cumsum(rnorm(length(t1))),k=10)
VMC<-VCM(cont=t1,theta=0.99)
sig<-runif(ncol(VMC))
VMC<-diag(sig)%*%VMC%*%diag(sig)
DTA<-mvrnorm(100,rep(0,ncol(VMC)),VMC)  

DTA<-sweep(DTA,2,CVC)
DTA<-apply(DTA,2,runmean,k=5)
matplot(t(DTA),type="l",col=1,lty=1)
user603
fonte
1
Você não pode apenas simular dados cuja média é uma função suave conhecida e adicionar ruído aleatório? Por exemplo,x=seq(0,2*pi,length=1000); plot(sin(x)+rnorm(1000)/10,type="l");
Macro
@ Macro: nop, se você aplicar mais zoom, verá que as funções geradas por ele não são suaves. Compare-os com algumas das curvas desses slides: bscb.cornell.edu/~hooker/FDA2007/Lecture1.pdf . Um spline suavizado do seu x poderia fazer o truque, mas estou procurando uma maneira direta de gerar os dados.
user603
sempre que você incluir ruído (que é uma parte necessária de qualquer modelo estocástico), os dados brutos serão, inerentemente, não suaves. O ajuste de spline a que você está se referindo está assumindo que o sinal é suave - não os dados reais observados (que é uma combinação de sinal e ruído).
Macro
@Macro: compare seus processos simulados para aqueles na página 16 deste documento: inference.phy.cam.ac.uk/mackay/gpB.pdf
user603
1
use polinômios de ordem superior. Um polinômio de 20º com coeficientes aleatórios (com a distribuição correta) pode mudar de direção (sem problemas) bastante. Se você encontrou uma resposta para sua pergunta, talvez possa publicá-la como resposta?
Macro

Respostas:

8

Veja como simular realizações de um Processo Gaussiano (GP). A suavidade das realizações depende das propriedades analíticas da função de covariância do GP. Este livro on-line possui muitas informações: http://unertaty.stat.cmu.edu/

Este vídeo fornece uma boa introdução aos GPs: http://videolectures.net/gpip06_mackay_gpb/

PS Em relação ao seu comentário, este código pode lhe dar um começo.

library(MASS)
C <- function(x, y) 0.01 * exp(-10000 * (x - y)^2) # covariance function
M <- function(x) sin(x) # mean function
t <- seq(0, 1, by = 0.01) # will sample the GP at these points
k <- length(t)
m <- M(t)
S <- matrix(nrow = k, ncol = k)
for (i in 1:k) for (j in 1:k) S[i, j] = C(t[i], t[j])
z <- mvrnorm(1, m, S)
plot(t, z)
zen
fonte
Você tem um link que aborda a questão de como simular realizações de um processo gaussiano, especificamente? Isso não é abordado no livro (olhando para o índice).
user603
A simulação de um GP é feita através das distribuições dimensionais finitas. Basicamente, você escolhe quantos pontos do domínio deseja e, a partir da função média e covariância do GP, obtém um normal multivariado. A amostragem dessa normal multivariada fornece o valor das realizações do GP nos pontos escolhidos. Como eu disse, esses valores se aproximam de uma função suave, desde que a função de covariância do GP satisfaça as condições analíticas necessárias. Uma função de covariância exponencial quadrática (com um termo "instável") é um bom começo.
Zen
4

{xEu,yEu}

require("MASS")
calcSigma<-function(X1,X2,l=1){
    Sigma<-matrix(rep(0,length(X1)*length(X2)),nrow=length(X1))
    for(i in 1:nrow(Sigma)){
        for (j in 1:ncol(Sigma)) Sigma[i,j]<-exp(-1/2*(abs(X1[i]-X2[j])/l)^2)
    }
    return(Sigma)
}
# The standard deviation of the noise
n.samples<-50
n.draws<-50
x.star<-seq(-5,5,len=n.draws)
nval<-3
f<-data.frame(x=seq(-5,5,l=nval),y=rnorm(nval,0,10))
sigma.n<-0.2
# Recalculate the mean and covariance functions
k.xx<-calcSigma(f$x,f$x)
k.xxs<-calcSigma(f$x,x.star)
k.xsx<-calcSigma(x.star,f$x)
k.xsxs<-calcSigma(x.star,x.star)
f.bar.star<-k.xsx%*%solve(k.xx+sigma.n^2*diag(1,ncol(k.xx)))%*%f$y
cov.f.star<-k.xsxs-k.xsx%*%solve(k.xx+sigma.n^2*diag(1,ncol(k.xx)))%*%k.xxs
values<-matrix(rep(0,length(x.star)*n.samples),ncol=n.samples)
for (i in 1:n.samples)  values[,i]<-mvrnorm(1,f.bar.star,cov.f.star)
values<-cbind(x=x.star,as.data.frame(values))
matplot(x=values[,1],y=values[,-1],lty=1,type="l",col="black")
lines(x.star,f.bar.star,col="red",lwd=2)

Um julgamento.  Funções suaves

user603
fonte
Este parece ser bom!
Zen