Jak symulować dane funkcjonalne?

12

Próbuję przetestować różne podejścia do analizy danych funkcjonalnych. Idealnie chciałbym przetestować zestaw podejść, jakie mam na symulowanych danych funkcjonalnych. Próbowałem wygenerować symulowany FD przy użyciu podejścia opartego na sumującym hałasie Gaussa (kod poniżej), ale uzyskane krzywe wyglądają o wiele za mocno w porównaniu do rzeczywistych .

Zastanawiałem się, czy ktoś ma wskaźnik do funkcji / pomysłów, aby wygenerować bardziej realistycznie wyglądające symulowane dane funkcjonalne. W szczególności powinny być gładkie. Jestem zupełnie nowy w tej dziedzinie, więc wszelkie porady są mile widziane.

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)
użytkownik603
źródło
1
Czy nie można po prostu symulować danych, których średnia jest znaną płynną funkcją i dodać losowy szum? Na przykładx=seq(0,2*pi,length=1000); plot(sin(x)+rnorm(1000)/10,type="l");
Macro
@Macro: nop, jeśli powiększysz wykres, zobaczysz, że generowane przez niego funkcje nie są płynne. Porównaj je z niektórymi krzywymi na tych slajdach: bscb.cornell.edu/~hooker/FDA2007/Lecture1.pdf . Wygładzony splajn twoich x może załatwić sprawę, ale szukam bezpośredniego sposobu na wygenerowanie danych.
user603
za każdym razem, gdy dołączasz szum (który jest niezbędną częścią każdego modelu stochastycznego), surowe dane z natury będą nierównomierne. Dopasowane do splajnu dopasowanie, o którym mówisz, zakłada, że sygnał jest płynny - a nie rzeczywiste obserwowane dane (które są kombinacją sygnału i szumu).
Makro
@Macro: porównaj symulowane procesy z tymi na stronie 16 tego dokumentu: inference.phy.cam.ac.uk/mackay/gpB.pdf
user603
1
używaj wielomianów wyższego rzędu. Wielomian dwudziestego stopnia z losowymi współczynnikami (z odpowiednim rozkładem) może całkiem (płynnie) zmieniać kierunki. Jeśli znalazłeś odpowiedź na swoje pytanie, być może możesz opublikować ją jako odpowiedź?
Makro,

Odpowiedzi:

8

Zobacz, jak symulować realizację procesu Gaussa (GP). Gładkość realizacji zależy od analitycznych właściwości funkcji kowariancji GP. Ta książka online zawiera wiele informacji: http://uncertainty.stat.cmu.edu/

Ten film stanowi miłe wprowadzenie do lekarzy rodzinnych: http://videolectures.net/gpip06_mackay_gpb/

PS Jeśli chodzi o komentarz, ten kod może dać ci początek.

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
źródło
Czy masz link, który odnosi się do pytania, jak konkretnie symulować realizację procesu Gaussa? Nie jest to omówione w książce (patrząc na indeks).
user603
Symulacja GP odbywa się poprzez skończone rozkłady wymiarowe. Zasadniczo wybierasz tyle punktów domeny, ile chcesz, a na podstawie średniej i funkcji kowariancji lekarza ogólnego uzyskujesz normalną wielowymiarową. Próbkowanie z tej wielowymiarowej normalnej daje wartość realizacji GP w wybranych punktach. Jak powiedziałem, wartości te są zbliżone do funkcji gładkiej, o ile funkcja kowariancji GP spełnia niezbędne warunki analityczne. Kwadratyczna funkcja kowariancji wykładniczej (z terminem „fluktuacji”) to dobry początek.
Zen,
4

{xi,yi}

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)

Próba.  Gładkie funkcje

użytkownik603
źródło
To wygląda dobrze!
Zen