# unlink("", recursive = TRUE)
#' Releasing Differential Private RKHS smoothing mean of a simulated dataset
#'
#' This function create a DP RKHS smoothing mean from an existing data set
#' with known eigenvalues and eigenvectors
#' @param alpha,beta Privacy parameters, real numbers
#' @param N real number, number of observations
#' @param n real number, number of grid points
#' @param e.val.x real vector n*1, eigenvalues
#' @param e.vec.x real valued matrix n*N, eigenvectors
#' @param tau range of the uniform distribution in KL expansion
#' @param phi real number, penalty parameter
#' @param mu real vector n*1, initial mean vector
#' @param e.val.z real valued matrix n*N, eigenvectors of noise
#' @param pow smoothing parameter, e.val.x_{i}=i^{-pow}
#' @param ro range parameter in kernel, real number
#' @return f.tilda: DP RKHS smoothing mean, n*1 real valued vector
#' @return delta: the coefficient of the noise, real number
#' @return f: RKHS smoothing mean, n*1 real valued vector
#' @return X: original data generated by the selected noise covariance operator C by KL epansion
#' when e.val.x_{i}=i^{-\code{pow}} and e.vec.x=e.vec.z ...
DP.Sim=function(alpha=1,beta=0.1,kernel="Gau",phi=0.01,ro=0.2,n=100,N=25,
tau=0.4,case=2,pow=4,mu=rep(0,n)){
grid=seq(0,1,length=n)
if(kernel=="Exp"){
C=function(t,s,ro){ # kernel of covariance operator of Exponential Process
return(exp(-abs(t-s)/ro))
}
}
if(kernel=="M3/2"){ # kernel of covariance operator of Matern Process nu=3/2
C=function(t,s,ro){
return((1+sqrt(3)*abs(t-s)/ro)*exp(-sqrt(3)*abs(t-s)/ro))
}
}
if(kernel=="M5/2"){ # kernel of covariance operator of Matern Process nu=5/2
C=function(t,s,ro){
return((1+(sqrt(5)*abs(t-s)/ro)+(5*(abs(t-s)^2)/(3*ro^2)))*exp(-sqrt(5)*abs(t-s)/ro))
}
}
if(kernel=="Gau"){ # kernel of covariance operator of Gaussian Process
C=function(t,s,ro){
return(exp(-(abs(t-s)^2)/ro))
}
}
Sig=matrix(nrow=n,ncol=n) # covariance matrix in the grid [0,1]
for(i in 1:n){ # Sig_{i,j}=C((i-1)/(n-1) , (j-1)/(n-1))
Sig[i,] = C(grid[i],grid,ro=ro)
}
####
gamma=eigen(Sig)$values
e.val.z=gamma[1:n]/n
e.vec.z=eigen(Sig)$vectors[,1:n]
m=length(which(e.val.z>0))
## generating Gaussian proccess Z based on KL-expansion
Z=matrix(0,nrow = n,ncol = 1)
for(i in 1:m){
Z=Z+sqrt(e.val.z[i])*rnorm(1)*e.vec.z[,i]
}
### generating f_{D} = muhat
if(case==0){
e.val.x=(e.val.z^2)
e.vec.x=e.vec.z
x=g.X(N=N,n=n,e.val.x=e.val.x,e.vec.x=e.vec.x,tau=tau,mu=mu)
f=rowMeans(x = x)
}
if(case==1){
B=seq(from = 1,to = length(grid))
e.val.x=(e.val.z*(B^-pow))
e.vec.x=e.vec.z
x=g.X(N=N,n=n,e.val.x=e.val.x,e.vec.x=e.vec.x,tau=tau,mu=mu)
f=rowMeans(x = x)
}
if(case==2){
B=seq(from = 1,to = length(grid))
e.val.x=as.vector(B^(-pow))
e.vec.x=e.vec.z
x=g.X(N=N,n=n,e.val.x=e.val.x,e.vec.x=e.vec.x,tau=tau,mu=mu,m=m) # it is needed for plots part
Y=g.X.RKHS(N=N,n=n,e.val.x=e.val.x,e.vec.x=e.vec.x,tau=tau,
phi=phi,mu=mu,m=m,e.val.z=e.val.z)
f=rowMeans(x = Y)
}
### calculating delta
if(case==0){
Delta2=(4*(tau^2)/N^2)*(sum(e.val.z))
delta=sqrt(((2*log(2/beta))/(alpha^2))*(Delta2))
}
if(case==1){
Delta2=(4*(tau^2)/N^2)*(sum(B^(-pow)))
delta=sqrt(((2*log(2/beta))/(alpha^2))*(Delta2))
}
if(case==2){
Delta2=(4*(tau^2)/N^2)*(sum(e.val.z*e.val.x/(e.val.z+phi)^2)) # phi should not be zero
delta=sqrt((2*log(2/beta)/(alpha^2))*Delta2)
}
### generating f.tilda
f.tilda=f+delta*Z
out=list("f.tilda"=f.tilda,"delta"=delta,"f"=f,"grid"=grid,"X"=x)
###
return(out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.