我最近遇到了高斯过程模型，并且碰巧认为它们可能是我在实验室中一直在研究的问题的解决方案。我有一个关于 Cross Validated 的开放且相关的问题，但我想将我的建模/数学问题与我的编程问题分开。因此，这是第二个相关的帖子。如果更多地了解我的问题的背景会有所帮助，尽管这里是我打开的CV 问题的链接。
这是我的 stan 代码，对应于我的 CV 帖子中提供的更新协方差函数：
functions{
//covariance function for main portion of the model
matrix main_GP(
int Nx,
vector x,
int Ny,
vector y,
real alpha1,
real alpha2,
real alpha3,
real rho1,
real rho2,
real rho3,
real rho4,
real rho5,
real HR_f,
real R_f){
matrix[Nx, Ny] K1;
matrix[Nx, Ny] K2;
matrix[Nx, Ny] K3;
matrix[Nx, Ny] Sigma;
//specifying random Gaussian process that governs covariance matrix
for(i in 1:Nx){
for (j in 1:Ny){
K1[i,j] = alpha1*exp(-square(x[i]-y[j])/2/square(rho1));
}
}
//specifying random Gaussian process incorporates heart rate
for(i in 1:Nx){
for(j in 1:Ny){
K2[i, j] = alpha2*exp(-2*square(sin(pi()*fabs(x[i]-y[j])*HR_f))/square(rho2))*
exp(-square(x[i]-y[j])/2/square(rho3));
}
}
//specifying random Gaussian process incorporates heart rate as a function of respiration
for(i in 1:Nx){
for(j in 1:Ny){
K3[i, j] = alpha3*exp(-2*square(sin(pi()*fabs(x[i]-y[j])*HR_f))/square(rho4))*
exp(-2*square(sin(pi()*fabs(x[i]-y[j])*R_f))/square(rho5));
}
}
Sigma = K1+K2+K3;
return Sigma;
}
//function for posterior calculations
vector post_pred_rng(
real a1,
real a2,
real a3,
real r1,
real r2,
real r3,
real r4,
real r5,
real HR,
real R,
real sn,
int No,
vector xo,
int Np,
vector xp,
vector yobs){
matrix[No,No] Ko;
matrix[Np,Np] Kp;
matrix[No,Np] Kop;
matrix[Np,No] Ko_inv_t;
vector[Np] mu_p;
matrix[Np,Np] Tau;
matrix[Np,Np] L2;
vector[Np] yp;
//--------------------------------------------------------------------
//Kernel Multiple GPs for observed data
Ko = main_GP(No, xo, No, xo, a1, a2, a3, r1, r2, r3, r4, r5, HR, R);
Ko = Ko + diag_matrix(rep_vector(1, No))*sn;
//--------------------------------------------------------------------
//kernel for predicted data
Kp = main_GP(Np, xp, Np, xp, a1, a2, a3, r1, r2, r3, r4, r5, HR, R);
Kp = Kp + diag_matrix(rep_vector(1, Np))*sn;
//--------------------------------------------------------------------
//kernel for observed and predicted cross
Kop = main_GP(No, xo, Np, xp, a1, a2, a3, r1, r2, r3, r4, r5, HR, R);
//--------------------------------------------------------------------
//Algorithm 2.1 of Rassmussen and Williams...
Ko_inv_t = Kop'/Ko;
mu_p = Ko_inv_t*yobs;
Tau=Kp-Ko_inv_t*Kop;
L2 = cholesky_decompose(Tau);
yp = mu_p + L2*rep_vector(normal_rng(0,1), Np);
return yp;
}
}
data {
int<lower=1> N1;
int<lower=1> N2;
vector[N1] X;
vector[N1] Y;
vector[N2] Xp;
real<lower=0> mu_HR;
real<lower=0> mu_R;
}
transformed data {
vector[N1] mu;
for(n in 1:N1) mu[n] = 0;
}
parameters {
real loga1;
real loga2;
real loga3;
real logr1;
real logr2;
real logr3;
real logr4;
real logr5;
real<lower=.5, upper=3.5> HR;
real<lower=.05, upper=.75> R;
real logsigma_sq;
}
transformed parameters {
real a1 = exp(loga1);
real a2 = exp(loga2);
real a3 = exp(loga3);
real r1 = exp(logr1);
real r2 = exp(logr2);
real r3 = exp(logr3);
real r4 = exp(logr4);
real r5 = exp(logr5);
real sigma_sq = exp(logsigma_sq);
}
model{
matrix[N1,N1] Sigma;
matrix[N1,N1] L_S;
//using GP function from above
Sigma = main_GP(N1, X, N1, X, a1, a2, a3, r1, r2, r3, r4, r5, HR, R);
Sigma = Sigma + diag_matrix(rep_vector(1, N1))*sigma_sq;
L_S = cholesky_decompose(Sigma);
Y ~ multi_normal_cholesky(mu, L_S);
//priors for parameters
loga1 ~ student_t(3,0,1);
loga2 ~ student_t(3,0,1);
loga3 ~ student_t(3,0,1);
logr1 ~ student_t(3,0,1);
logr2 ~ student_t(3,0,1);
logr3 ~ student_t(3,0,1);
logr4 ~ student_t(3,0,1);
logr5 ~ student_t(3,0,1);
logsigma_sq ~ student_t(3,0,1);
HR ~ normal(mu_HR,.25);
R ~ normal(mu_R, .03);
}
generated quantities {
vector[N2] Ypred;
Ypred = post_pred_rng(a1, a2, a3, r1, r2, r3, r4, r5, HR, R, sigma_sq, N1, X, N2, Xp, Y);
}
我已经修改了我的内核中包含的参数的先验，一些参数化有点快（在某些情况下快两倍，但即使在这些情况下仍然可以产生相对较慢的链）。
我正在尝试使用受污染部分前后 15 秒的数据（以 3.33 Hz 采样，因此总共 100 个数据点）来预测价值 3.5 秒的数据（以 10 Hz 采样 - 所以 35 个数据点）的值。
R中指定的模型如下：
fit.pred2 <- stan(file = 'Fast_GP6_all.stan',
data = dat,
warmup = 1000,
iter = 1500,
refresh=5,
chains = 3,
pars= pars.to.monitor
)
老实说，我不知道我是否需要那么多热身迭代。我想部分缓慢的估计是相当无信息的先验结果（心率和呼吸除外，HR因为R它们在健康成年人的休息时具有相当广为人知的范围）。
非常欢迎任何建议来加快我的程序的运行时间。
谢谢。