我正在编写一些代码来从贝叶斯的角度运行高斯过程，并且我想在循环中估计我的参数，然后还要在循环中获得克里金估计器。我遇到的问题是我必须在 R 中运行一个双 for 循环，这可能非常慢。我试图弄清楚如何加快计算速度（尝试将克里金估计器编写为函数，然后将其传递给循环内的应用函数），但它仍然非常缓慢。希望有人能够就如何矢量化我的问题或任何其他可能加快代码速度的技巧提供一些指导。这是我写的代码：
#Example Data and function
x = sort(runif(4,0,1))
y = exp(-1.4*x)*cos(7*pi*x/2)
#Distance matrix
tau = as.matrix(dist(x,upper=T,diag=T))
#Correlation Matrix
corrR = function(psi,tau){
R = exp(-tau^(1.99)*psi)
return(R)
}
#Full conditional for psi
psi.cond = function(psi,r,sig,D,beta,Y){
d = .01
e = .01
ans = det(r)^(-.5)*exp(-.5/sig*t(Y-D*beta)%*%solve(r)%*%(Y-D*beta))*psi^-(d+1)*exp(-e/psi)
#ans = -.5*log(det(r))-.5/sig*t(Y-D*beta)%*%solve(r)%*%(Y-D*beta)-(d+1)*log(psi)-e/psi
return(as.real(ans))
}
#Kriging estimator
krig = function(x.new,x,beta,psi,tau,Y,D){
S = length(x.new)
pred = rep(NA,S)
for(j in 1:S){
R = corrR(psi,tau)
r = as.matrix(dist(c(x,x.new[j]),upper=T,diag=T))[(length(x)+1):(length(x)+length(x.new[j])),1:length(x)]
r = t(as.matrix(exp(-psi*r^(1.99)),nrow=nrow(r),ncol=length(x)))
pred[j] = as.real(beta+r%*%solve(R)%*%(Y-D*beta))
}
return(pred)
}
D = rep(1,length(x))
Y = as.matrix(y)
m = length(x)
a = 2
b = 1
#Number of MCMC iterations = B
B = 50000
beta = c(1,rep(NA,B))
sigma = rep(NA,B)
psi = c(120,rep(NA,B-1))
#Number of predicted points = S
S = 100
yhat = matrix(NA,nrow=B,ncol=S)
x.new = as.matrix(seq(0,1,len=S))
for(i in 1:B){
R = corrR(psi[i],tau)
bhat = as.real(solve(t(D)%*%solve(R)%*%D)%*%t(D)%*%solve(R)%*%Y)
sigma[i] = 1/rgamma(1,(m+2*a)/2,(as.real(t(Y-D*beta[i])%*%solve(R)%*%(Y-D*beta[i]))+2*b)/2)
beta[i+1] = rnorm(1,bhat,t(D)%*%solve(sigma[i]*R)%*%D)
log.xi = rnorm(1,log(psi[i]),.1)
xi = exp(log.xi)
u = runif(1)
R.xi = corrR(xi,tau)
R.psi = corrR(psi[i],tau)
temp = (psi.cond(xi,R.xi,sigma[i],D,beta[i],Y)*(1/psi[i]))/(psi.cond(psi[i],R.psi,sigma[i],D,beta[i],Y)*(1/xi))
alpha = min(1,temp)
if(u <= alpha){
psi[i+1] = xi
}else{
psi[i+1] = psi[i]
}
yhat[i,] = apply(x.new,1,krig,x=x,beta=beta[i+1],psi=psi[i+1],tau,Y,D)
}