R/DP.Sim.R

In sxz155/PFDA:

# 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)
}