#----------------------------------------------------------------------------------------#
#---------------------------- CONTINUOUS OBSERVATIONS -----------------------------------#
#----------------------------------------------------------------------------------------#
# smc2 version when predictive density is intractable
hincrement_continuous_smc2_ = function(t,model,observationt,thetas,Wtheta,PFs) {
Ntheta = ncol(thetas)
hincrement = 0
d1log = list()
d2log = list()
for (k in 1:model$dimY) {
Ek_theta = vector("numeric",Ntheta)
Fk_theta = vector("numeric",Ntheta)
for (m in 1:Ntheta){
Ek_theta[m] = 0
Fk_theta[m] = 0
if (Wtheta[m] == 0) {
next
}
PF = PFs[[m]]
X = PF$X
WX = PF$xnormW
if (is.null(dim(X))){
X = matrix(X,nrow = model$dimX)
}
# compute the derivatives for all k and all X at once during the first loop
if (k==1){
derivatives = model$derivativelogdobs(observationt,X,t,thetas[,m])
d1log[[m]] = derivatives$jacobian
d2log[[m]] = derivatives$hessiandiag
}
# we automatically discard particles with weight 0 to avoid artificial apparition of NaN
# (which could happen when performing 0*Inf)
nonzeroWX = (WX>0)
Ek_theta[m] = sum(WX[nonzeroWX]*(d2log[[m]][nonzeroWX,k] + (d1log[[m]][nonzeroWX,k])^2))
Fk_theta[m] = sum(WX[nonzeroWX]*d1log[[m]][nonzeroWX,k])
}
# we automatically discard particles with weight 0 to avoid artificial apparition of NaN
# (which could happen when performing 0*Inf)
nonzeroWtheta = (Wtheta>0)
Ek = sum(Wtheta[nonzeroWtheta]*Ek_theta[nonzeroWtheta])
Fk = (sum(Wtheta[nonzeroWtheta]*Fk_theta[nonzeroWtheta]))^2
hincrement = hincrement + 2*Ek - Fk
}
return (hincrement)
}
#-------------------------------------------------------------------------------------------
# smc version when predictive density is available
hincrement_continuous_smc_ = function(t,model,observations,thetas,Wtheta,byproducts) {
Ntheta = ncol(thetas)
hincrement = 0
d1log = list()
d2log = list()
for (k in 1:model$dimY) {
Ek_theta = array(0,dim = Ntheta)
Fk_theta = array(0,dim = Ntheta)
for (m in 1:Ntheta){
if (Wtheta[m] == 0) {
next
}
# compute the derivatives for all k at once during the first loop
if (k==1){
if (!is.null(byproducts)){
derivatives = model$derivativelogdpredictive(observations,t,thetas[,m],byproducts[[m]])
} else {
derivatives = model$derivativelogdpredictive(observations,t,thetas[,m])
}
d1log[[m]] = derivatives$jacobian
d2log[[m]] = derivatives$hessiandiag
}
Ek_theta[m] = d2log[[m]][k] + (d1log[[m]][k])^2
Fk_theta[m] = d1log[[m]][k]
}
# we automatically discard particles with weight 0 to avoid artificial apparition of NaN
# (which could happen when performing 0*Inf)
nonzeroWtheta = (Wtheta>0)
Ek = sum(Wtheta[nonzeroWtheta]*Ek_theta[nonzeroWtheta])
Fk = (sum(Wtheta[nonzeroWtheta]*Fk_theta[nonzeroWtheta]))^2
hincrement = hincrement + 2*Ek - Fk
}
return (hincrement)
}
#----------------------------------------------------------------------------------------#
#------------------------------ DISCRETE OBSERVATIONS -----------------------------------#
#----------------------------------------------------------------------------------------#
# This function computes the approximation pt_hat(y) for a given set of particles
phat_smc2 = function(t,model,y,thetas,thetanormw,Ntheta,Xpred,xprednormw) {
py = 0
for (m in 1:Ntheta){
if (thetanormw[m]==0){
next
}
py = py + thetanormw[m]*sum(xprednormw[,m]*model$dobs(y,Xpred[,,m],t,thetas[,m],log = FALSE))
}
return (py)
}
# This function computes the partial score term Hdk (forward difference version)
# a,b are vectors of componentwise lower and upper bounds of the observations
Hdk_forward = function(k,a,b,d,y,py_minusek,py,py_plusek) {
if (y[k]==b[k]) {
return (-2*(py-py_minusek)/py_minusek)
}
else {
if (y[k]==a[k]) {
return (2*(py_plusek-py)/py + ((py_plusek-py)/py)^2)
}
else {
return (2*((py_plusek-py)/py-(py-py_minusek)/py_minusek) + ((py_plusek-py)/py)^2)
}
}
}
# This function computes the partial score term Hdk (central difference version)
# a,b are vectors of componentwise lower and upper bounds of the observations
Hdk_central = function(k,a,b,d,y,py_minus2ek,py_minusek,py,py_plusek,py_plus2ek) {
if (y[k]==b[k]) {
return (-(py-py_minus2ek)/(2*py_minusek))
}
if (y[k]==b[k]-1) {
return (-(py-py_minus2ek)/(2*py_minusek) + ((py_plusek-py_minusek)/(2*py))^2)
}
if (y[k]==a[k]) {
return ((py_plus2ek-py)/(2*py_plusek))
}
if (y[k]==a[k]+1) {
return ((py_plus2ek-py)/(2*py_plusek) + ((py_plusek-py_minusek)/(2*py))^2)
}
if (y[k]!=b[k] && y[k]!=b[k]-1 && y[k]!=a[k] && y[k]!=a[k]+1){
return ((py_plus2ek-py)/(2*py_plusek) - (py-py_minus2ek)/(2*py_minusek) + ((py_plusek-py_minusek)/(2*py))^2)
}
}
# This function computes the partial score term Hdk and is a wrapper of Hdk_forward and Hdk_central
# a,b are vectors of componentwise lower and upper bounds of the observations
Hdk = function(k,a,b,d,y,py_minusek,py,py_plusek,diff_type = "central",py_minus2ek = NULL,py_plus2ek = NULL) {
if (diff_type == "central"){
return (Hdk_central(k,a,b,d,y,py_minus2ek,py_minusek,py,py_plusek,py_plus2ek))
}
if (diff_type == "forward"){
return (Hdk_forward(k,a,b,d,y,py_minusek,py,py_plusek))
}
}
# This function computes the partial score term Hd
# a,b are vectors of componentwise lower and upper bounds of the observations
# d is the dimension of y
Hd_smc2 = function(t,model,yt,thetas,thetanormw,Xpred,xprednormw, diff_type = "central") {
Ntheta = ncol(thetas)
a = model$lower
b = model$upper
d = model$dimY
result = 0
for (k in 1:d) {
ek = rep(0,d)
ek[k] = 1
py = phat_smc2(t,model,yt,thetas,thetanormw,Ntheta,Xpred,xprednormw)
py_minusek = phat_smc2(t,model,yt-ek,thetas,thetanormw,Ntheta,Xpred,xprednormw)
py_plusek = phat_smc2(t,model,yt+ek,thetas,thetanormw,Ntheta,Xpred,xprednormw)
if (diff_type == "central") {
py_minus2ek = phat_smc2(t,model,yt-2*ek,thetas,thetanormw,Ntheta,Xpred,xprednormw)
py_plus2ek = phat_smc2(t,model,yt+2*ek,thetas,thetanormw,Ntheta,Xpred,xprednormw)
} else {
py_minus2ek = NULL
py_plus2ek = NULL
}
result = result + Hdk(k,a,b,d,yt,py_minusek,py,py_plusek, diff_type, py_minus2ek, py_plus2ek)
}
return (result)
}
#-------------------------------------------------------------------------------------------
# smc version when predictive density is available
# This function computes the approximation qt_hat(y) for a given set of particles
phat_smc = function(t,model,observations,thetas,thetanormw,Ntheta,byproducts) {
py = 0
if (!is.null(byproducts)){
for (m in 1:Ntheta){
if (thetanormw[m]==0){
next
}
py = py + thetanormw[m]*model$dpredictive(observations,t,thetas[,m],byproducts[[m]],log = FALSE)
}
} else {
for (m in 1:Ntheta){
if (thetanormw[m]==0){
next
}
py = py + thetanormw[m]*model$dpredictive(observations,t,thetas[,m],log = FALSE)
}
}
return (py)
}
# This function computes the partial score term Hd
# a,b are vectors of componentwise lower and upper bounds of the observations
# d is the dimension of y
Hd_smc = function(t,model,observations,thetas,thetanormw,byproducts, diff_type = "central") {
Ntheta = ncol(thetas)
a = model$lower
b = model$upper
d = model$dimY
result = 0
for (k in 1:d) {
ek = rep(0,d)
ek[k] = 1
obsminusek = observations
obsplusek = observations
obsminusek[,t] = obsminusek[,t]-ek
obsplusek[,t] = obsplusek[,t]+ek
if (diff_type == "central") {
obsminus2ek = observations
obsplus2ek = observations
obsminus2ek[,t] = obsminus2ek[,t]-2*ek
obsplus2ek[,t] = obsplus2ek[,t]+2*ek
} else {
obsminus2ek = NULL
obsplus2ek = NULL
}
py = phat_smc(t,model,observations,thetas,thetanormw,Ntheta,byproducts)
py_minusek = phat_smc(t,model,obsminusek,thetas,thetanormw,Ntheta,byproducts)
py_plusek = phat_smc(t,model,obsplusek,thetas,thetanormw,Ntheta,byproducts)
if (diff_type == "central") {
py_minus2ek = phat_smc(t,model,obsminus2ek,thetas,thetanormw,Ntheta,byproducts)
py_plus2ek = phat_smc(t,model,obsplus2ek,thetas,thetanormw,Ntheta,byproducts)
} else {
py_minus2ek = NULL
py_plus2ek = NULL
}
result = result + Hdk(k,a,b,d,observations[,t],py_minusek,py,py_plusek, diff_type, py_minus2ek, py_plus2ek)
}
return (result)
}
#-------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------
# Wrappers: the following function performs all the preprocessing needed to compute
# the Hyvarinen score (e.g. generate x-particles from the one-step predictive distribution,
# generate more theta-particles for variance reduction if needed, etc ...)
#-------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------
# Discrete case, SMC2: wrapper of Hd_smc2
hincrement_discrete_smc2 = function(thetas, normw, PFs, t, observations, model,
logtargetdensities, algorithmic_parameters) {
discrete_diff_type = algorithmic_parameters$discrete_diff_type
# If reduce_variance: temporarily increase the number of particles before computing
# the Hyvarinen score.
if (algorithmic_parameters$reduce_variance) {
# generate additional theta-particles (and their associated particle filters)
larger_pool = get_additional_particles_smc2(algorithmic_parameters$Nc, thetas, normw, PFs, t,
observations, model, logtargetdensities,
algorithmic_parameters, algorithmic_parameters$Ncx)
thetas_pool = larger_pool$thetas
PFs_pool = larger_pool$PFs
normw_pool = rep(1/ncol(thetas_pool), ncol(thetas_pool))
# Get the new number of particles
Ntheta_pool = ncol(thetas_pool)
Nx_pool = PFs_pool[[1]]$Nx
# generate x-particles targeting the one-step predictive distribution which are already
# contained in the PFs after initialization
if (t==1) {
Xpred = array(NA,dim = c(model$dimX,Nx_pool,Ntheta_pool)) # (need to reconstruct since size Nx might change)
XnormW_previous = matrix(1/Nx_pool, nrow = Nx_pool, ncol = Ntheta_pool)
for (itheta in 1:Ntheta_pool){
Xpred[,,itheta] = PFs_pool[[itheta]]$X
}
}
# Construct particles targeting the one-step-ahead predictive
if (t > 1){
Xpred = array(NA,dim = c(model$dimX,Nx_pool,Ntheta_pool))
XnormW_previous = matrix(NA, nrow = Nx_pool, ncol = Ntheta_pool)
# At this point in the code, PFs[[itheta]]$X contains the x-particles from the filtering
# distribution at time (t-1)
for (itheta in 1:Ntheta_pool){
if (is.null(dim(PFs_pool[[itheta]]$X))){
XnormW_previous[,itheta] = PFs_pool[[itheta]]$xnormW
Xpred[,,itheta] = model$rtransition(matrix(PFs_pool[[itheta]]$X,ncol = Nx_pool), t, thetas_pool[,itheta])
}
else{
XnormW_previous[,itheta] = PFs_pool[[itheta]]$xnormW
Xpred[,,itheta] = model$rtransition(PFs_pool[[itheta]]$X, t, thetas_pool[,itheta])
}
}
}
# compute incremental H score (with theta from time t-1, see formula in the paper)
return (Hd_smc2(t,model,observations[,t],thetas_pool,normw_pool,Xpred,XnormW_previous, discrete_diff_type))
} else {
Ntheta = ncol(thetas)
Nx = PFs[[1]]$Nx
# Construct particles targeting the one-step-ahead predictive (need to reconstruct since size Nx might change)
if (t==1) {
Xpred = array(NA,dim = c(model$dimX, Nx, Ntheta))
XnormW_previous = matrix(1/Nx, nrow = Nx, ncol = Ntheta) #matrix of normalized weights for X at previous step
for (itheta in 1:Ntheta){
Xpred[,,itheta] = PFs[[itheta]]$X
}
}
if (t > 1){
Nx = PFs[[1]]$Nx
Xpred = array(NA,dim = c(model$dimX,Nx,Ntheta)) # (need to reconstruct since size Nx might change)
XnormW_previous = matrix(NA, nrow = Nx, ncol = Ntheta)
# At this point in the code, PFs[[itheta]]$X contains the x-particles from the filtering
# distribution at time (t-1)
for (itheta in 1:Ntheta){
if (is.null(dim(PFs[[itheta]]$X))){
XnormW_previous[,itheta] = PFs[[itheta]]$xnormW
Xpred[,,itheta] = model$rtransition(matrix(PFs[[itheta]]$X,ncol = Nx), t, thetas[,itheta])
}
else{
XnormW_previous[,itheta] = PFs[[itheta]]$xnormW
Xpred[,,itheta] = model$rtransition(PFs[[itheta]]$X, t, thetas[,itheta])
}
}
}
# compute incremental H score (with theta from time t-1, see formula in the paper)
return (Hd_smc2(t,model,observations[,t],thetas,normw,Xpred,XnormW_previous, discrete_diff_type))
}
}
#-------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------
# Discrete case, SMC: wrapper of Hd_smc
hincrement_discrete_smc = function(thetas, normw, byproducts, t, observations, model,
logtargetdensities, algorithmic_parameters) {
discrete_diff_type = algorithmic_parameters$discrete_diff_type
# compute incremental H score (with theta from time t-1, see formula in the paper)
if (algorithmic_parameters$reduce_variance) {
# generate additional particles
larger_pool = get_additional_particles_smc(algorithmic_parameters$Nc, thetas, normw,
byproducts, t, observations, model,
logtargetdensities, algorithmic_parameters)
thetas_pool = larger_pool$thetas
byproducts_pool = larger_pool$byproducts
normw_pool = rep(1/ncol(thetas_pool), ncol(thetas_pool))
# compute Hyvarinen score with more particles
return (Hd_smc(t,model,observations,thetas_pool,normw_pool,byproducts_pool,diff_type = discrete_diff_type))
} else {
return (Hd_smc(t,model,observations,thetas,normw,byproducts,diff_type = discrete_diff_type))
}
}
#-------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------
# Continuous case, SMC2: wrapper of hincrement_continuous_smc2_
hincrement_continuous_smc2 = function(thetas, normw, PFs, t, observations, model,
logtargetdensities, algorithmic_parameters) {
# If reduce_variance: temporarily increase the number of particles before computing
# the Hyvarinen score.
if (algorithmic_parameters$reduce_variance) {
# generate additional theta-particles (and their associated particle filters)
larger_pool = get_additional_particles_smc2(algorithmic_parameters$Nc, thetas, normw, PFs, t,
observations, model, logtargetdensities,
algorithmic_parameters, algorithmic_parameters$Ncx)
thetas_pool = larger_pool$thetas
PFs_pool = larger_pool$PFs
normw_pool = rep(1/ncol(thetas_pool), ncol(thetas_pool))
return (hincrement_continuous_smc2_(t,model,observations[,t,drop=FALSE],thetas_pool,normw_pool,PFs_pool))
} else {
return (hincrement_continuous_smc2_(t,model,observations[,t,drop=FALSE],thetas,normw,PFs))
}
}
#-------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------
# Continuous case, SMC: wrapper of hincrement_continuous_smc_
hincrement_continuous_smc = function(thetas, normw, byproducts, t, observations, model,
logtargetdensities, algorithmic_parameters) {
# If reduce_variance: temporarily increase the number of particles before computing
# the Hyvarinen score.
if (algorithmic_parameters$reduce_variance) {
# generate additional particles
larger_pool = get_additional_particles_smc(algorithmic_parameters$Nc, thetas, normw,
byproducts, t, observations, model,
logtargetdensities, algorithmic_parameters)
thetas_pool = larger_pool$thetas
byproducts_pool = larger_pool$byproducts
normw_pool = rep(1/ncol(thetas_pool), ncol(thetas_pool))
# compute Hyvarinen score with more particles
return (hincrement_continuous_smc_(t, model, observations,thetas_pool,normw_pool,byproducts_pool))
} else {
return (hincrement_continuous_smc_(t, model, observations,thetas,normw,byproducts))
}
}
#-------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------
# Additional implementations using Kernel density estimators
#-------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------
# smc version when predictive density is available, using density estimators
# via local regression. WARNING: only implemented for univariate observations.
hincrementContinuous_smc_dde = function(t,model,observations,thetas,normw,byproducts,
logtargetdensities, algorithmic_parameters){
if (t > algorithmic_parameters$dde_options$nb_steps){
return (NA)
}
if (is.null(model$rpredictive)) {
print("Can't perform density estimation: no sampler 'rpredictive' provided")
return (NA)
}
if (model$dimY!=1) {
print("WARNING: density estimation has only been implemented for univariate observations")
return (NA)
} else {
Ny = algorithmic_parameters$dde_options$Ny
sigma20 = algorithmic_parameters$dde_options$sigma2_order0
sigma21 = algorithmic_parameters$dde_options$sigma2_order1
sigma22 = algorithmic_parameters$dde_options$sigma2_order2
if (t == 1){
thetas_sim = model$rprior(Ny)
Yt_sim = apply(thetas_sim,MARGIN = 2, function(theta)model$rpredictive(1,1,theta,NULL))
}
if (t >=2){
# generate additional particles until we have at least Ny particles
larger_pool = get_additional_particles_smc(Ny, thetas, normw, byproducts, t, observations, model,
logtargetdensities, algorithmic_parameters)
# Note that the particles thetas returned by get_additional_particles_smc are equally weighted
thetas_pool = larger_pool$thetas
# Draw one Yt for each particle theta
Yt_sim = apply(thetas_pool,MARGIN = 2, function(theta)model$rpredictive(1,t,theta,observations[,1:(t-1),drop=FALSE]))
}
d0_est = get_derivative_RBFlocal(Yt_sim,observations[,t],sigma20,order = 0)
d1_est = get_derivative_RBFlocal(Yt_sim,observations[,t],sigma21,order = 1)
d2_est = get_derivative_RBFlocal(Yt_sim,observations[,t],sigma22,order = 2)
# compute h score via local estimation
return ((2*d2_est/d0_est)-((d1_est/d0_est)^2))
}
}
# smc version when predictive density is available, using density estimators
# via local regression. WARNING: only implemented for univariate observations.
hincrementContinuous_smc2_dde = function(t,model,observations,thetas,normw,PFs,
logtargetdensities, algorithmic_parameters){
if (t > algorithmic_parameters$dde_options$nb_steps){
return (NA)
}
if (is.null(model$robs)) {
print("Can't perform density estimation: no sampler 'robs' provided")
return (NA)
}
if (model$dimY!=1) {
print("WARNING: density estimation has only been implemented for univariate observations")
return (NA)
} else {
Ny = algorithmic_parameters$dde_options$Ny
sigma20 = algorithmic_parameters$dde_options$sigma2_order0
sigma21 = algorithmic_parameters$dde_options$sigma2_order1
sigma22 = algorithmic_parameters$dde_options$sigma2_order2
if (t == 1){
thetas_pool = model$rprior(Ny)
Xt_sim = apply(thetas_pool,MARGIN = 2, function(theta)model$rinitial(theta,1))
Ntheta_larger = ncol(thetas_pool)
Yt_sim = rep(NA,Ntheta_larger)
for (i in 1:Ntheta_larger){
Yt_sim[i] = model$robs(Xt_sim[,i,drop=FALSE],1,thetas_pool[,i,drop=FALSE])
}
}
if (t >=2){
# generate additional particles until we have at least Ny particles
larger_pool = get_additional_particles_smc2(Ny, thetas, normw, PFs, t, observations, model,
logtargetdensities, algorithmic_parameters, NULL)
# Note that the particles thetas returned by get_additional_particles_smc are equally weighted
thetas_pool = larger_pool$thetas
PFs_pool = larger_pool$PFs
Ntheta_larger = ncol(thetas_pool)
Yt_sim = rep(NA,Ntheta_larger)
for (i in 1:Ntheta_larger){
Xprev = PFs_pool[[i]]$X[,sample(1:PFs_pool[[i]]$Nx,size = 1,prob = PFs_pool[[i]]$xnormW),drop=FALSE]
Xt_sim = model$rtransition(Xprev,t,thetas_pool[,i])
Yt_sim[i] = model$robs(Xt_sim,t,thetas_pool[,i])
}
}
d0_est = get_derivative_RBFlocal(Yt_sim,observations[,t],sigma20,order = 0)
d1_est = get_derivative_RBFlocal(Yt_sim,observations[,t],sigma21,order = 1)
d2_est = get_derivative_RBFlocal(Yt_sim,observations[,t],sigma22,order = 2)
# compute h score via local estimation
return ((2*d2_est/d0_est)-((d1_est/d0_est)^2))
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.