Как симулировать функциональные данные?

12

Я пытаюсь проверить различные подходы анализа функциональных данных. В идеале я хотел бы протестировать панель подходов, которые у меня есть, на смоделированных функциональных данных. Я попытался сгенерировать смоделированный FD, используя подход, основанный на суммировании гауссовских шумов (код ниже), но полученные кривые выглядят слишком грубыми по сравнению с реальными .

Мне было интересно, есть ли у кого-нибудь указатель на функции / идеи, чтобы генерировать более реалистично выглядящие моделируемые функциональные данные. В частности, они должны быть гладкими. Я совершенно новичок в этой области, поэтому любые советы приветствуются.

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
источник
1
Разве вы не можете просто смоделировать данные, среднее значение которых является известной гладкой функцией, и добавить случайный шум? Например,x=seq(0,2*pi,length=1000); plot(sin(x)+rnorm(1000)/10,type="l");
Макро
@Macro: нет, если вы увеличите график, вы увидите, что сгенерированные им функции не являются плавными. Сравните их с некоторыми кривыми на этих слайдах: bscb.cornell.edu/~hooker/FDA2007/Lecture1.pdf . Сглаженный сплайн ваших х может помочь, но я ищу прямой способ генерации данных.
user603
всякий раз, когда вы включаете шум (который является необходимой частью любой стохастической модели), необработанные данные по своей природе будут негладкими. Подход сплайна, на который вы ссылаетесь, предполагает, что сигнал ровный, а не фактические наблюдаемые данные (которые представляют собой комбинацию сигнала и шума).
Макро
@Macro: сравните ваши смоделированные процессы с теми, что приведены на странице 16 этого документа: inference.phy.cam.ac.uk/mackay/gpB.pdf
user603
1
использовать многочлены высшего порядка. Полином 20-й степени со случайными коэффициентами (с правильным распределением) может довольно сильно менять направления (плавно). Если вы нашли ответ на свой вопрос, возможно, вы можете опубликовать его как ответ?
Макро

Ответы:

8

Посмотрите, как моделировать реализации гауссовского процесса (GP). Гладкость реализаций зависит от аналитических свойств ковариационной функции ГП. В этой онлайн-книге много информации: http://unterminty.stat.cmu.edu/

Это видео дает хорошее представление о врачах общей практики: http://videolectures.net/gpip06_mackay_gpb/

PS Что касается вашего комментария, этот код может дать вам начало.

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
источник
У вас есть ссылка, которая решает вопрос о том, как, в частности, моделировать реализации гауссовского процесса? Об этом не говорится в книге (глядя на указатель).
user603
Моделирование ГП осуществляется через конечномерные распределения. По сути, вы выбираете столько точек домена, сколько хотите, и из средней и ковариационной функции GP вы получаете многовариантную нормаль. Выборка из этого многомерного нормального дает вам значение реализации ГП в выбранных точках. Как я уже сказал, эти значения аппроксимируют гладкую функцию, если ковариационная функция GP удовлетворяет необходимым аналитическим условиям. Квадратичная экспоненциальная ковариационная функция (с термином «джиттер») является хорошим началом.
Дзен
4

Хорошо, вот ответ, который я придумал (по сути, он взят здесь и здесь ). Идея состоит в том, чтобы спроецировать несколько случайных пар на сплайн-основу. Тогда мы уверены, что получим ничью от (гладкого) GP.{Икся,Yя}

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)

Испытание.  Гладкие функции

user603
источник
Это выглядит хорошо!
Дзен