R/RcppExports.R

In SunJinglan/StatComp21080: Density estimations using Lindsey's method and penalized Gaussian mixtures

# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

#' @title Lindsey's Method
#' @description A density estimation using Lindsey's method with Gibbs sampling
#' @param df the sample to be estimated
#' @param n the number of equal width bins
#' @param a the minimum of the range to estimate the density
#' @param b the maximum of the range to estimate the density
#' @param iter the number of Gibbs samples
#' @param burnin the number of between-sample random numbers
#' @return a list where \code{t} represents the possilble values and \code{lm} represents the density of that value 
#' @examples
#' \dontrun{
#' rnR <- gibbsR(100,10)
#' par(mfrow=c(2,1));
#' plot(rnR[,1],type='l')
#' plot(rnR[,2],type='l')
#' }
#' @import dplyr
#' @import SuppDists
#' @import MASS
#' @export
LM<-function(df,n=20,a=min(df),b=max(df),iter=5000,burnin=1000){
  #df：估计对象
  #n：分段数
  #N,thin：吉布斯采样次数
  # 模型适应 
  nn<-length(df)
  k<-n
  abk<-seq(a,b,length.out = (k+1))
  t<-(abk[-1]+abk[-(k+1)])/2
  N<-table(cut(df,abk))
  Y<-matrix(sqrt(k/n*(N+1/4)))
  X<-matrix(c(rep(1,n),t),ncol=2)
  Omega<-matrix(0, ncol = n,nrow = n)
  for (i in 1:n) {
    for (j in 1:n) {
      if (t[i]==t[j]) Omega[i,j]<-1
      if (t[i]<t[j]) Omega[i,j]<-t[i]^2*(t[j]-t[i]/3)/2
      if (t[i]>t[j]) Omega[i,j]<-t[j]^2*(t[i]-t[j]/3)/2
    }
  }
  D<-diag(eigen(Omega)$values)
  Q<-eigen(Omega)$vectors
  m<-sum(D>0)
  D<-D[,1:m]
  Z<-Q%*%sqrt(D)
  XZ<-cbind(X,Z)
  
  # Gibbs采样
  c.tau<-10^5
  c.sigma<-10^3
  sigma.beta<-10^8
  tau<-rep(1,iter)
  sigma<-rep(1,iter)
  beta.u<-matrix(1,ncol=2+m,nrow=iter)
  for (i in 2:iter) {
    while (TRUE) {
      tau[i]<-1/rgamma(1,shape=m/2-1,rate=beta.u[i-1,3:(2+m)]%*%beta.u[i-1,3:(2+m)]/2)
      if(tau[i]<=c.tau) break
    }
    A<-diag(c(rep(sigma.beta,2),rep(tau[i],m)))
    while (TRUE) {
      sigma[i]<-1/rgamma(1,shape=n/2-1,rate=t(Y-XZ%*%beta.u[i-1,])%*%(Y-XZ%*%beta.u[i-1,])/2)
      if(sigma[i]<=c.sigma) break
    }
    beta.u[i,]<-mvrnorm(mu=solve(t(XZ)%*%XZ+sigma[i]*solve(A))%*%t(XZ)%*%Y,
                        Sigma=sigma[i]*solve(t(XZ)%*%XZ+sigma[i]*solve(A)))
  }
  b<-burnin+1
  tau.hat<-mean(tau[b:iter])
  sigma.hat<-mean(sigma[b:iter])
  beta.hat<-colMeans(beta.u[b:iter,1:2])
  u.hat<-colMeans(beta.u[b:iter,3:(2+m)])
  
  # 模型代入
  r.hat<-beta.hat[1]+beta.hat[2]*t+Z%*%u.hat
  r.p<-apply(matrix(c(rep(0,n),r.hat),ncol=2),1, max)
  I<-sum((t[2:n]-t[1:(n-1)])*(r.p[2:n]^2+r.p[1:(n-1)]^2)/2)
  lm<-r.p^2/I
  return(list(t=t,lm=lm))
}


#' @title Penalized Gaussian Mixtures
#' @description A density estimation using Penalized Gaussian Mixtures method with Gibbs sampling and Hamiltonian Monte Carlo sampling
#' @param x the sample to be estimated
#' @param K the number of Gaussian distributions to be cumulated
#' @param a0 the minimum of the range to estimate the density
#' @param b0 the maximum of the range to estimate the density
#' @param iter the number of Gibbs samples
#' @param burnin the number of between-sample random numbers
#' @return a list where \code{t} represents the possilble values and \code{pgm} represents the density of that value 
#' @examples
#' \dontrun{
#' rnR <- gibbsR(100,10)
#' par(mfrow=c(2,1));
#' plot(rnR[,1],type='l')
#' plot(rnR[,2],type='l')
#' }
#' @import mixtools 
#' @export
PGM<-function(x,K=35,a0=min(x),b0=max(x),iter=5000,burnin=1000){#x:样本,a.b:上下限,k:分段数
  #模型适应
  n<-length(x)
  abk<-seq(a0,b0,length.out = (K+1))
  mu<-(abk[-1]+abk[-(K+1)])/2
  delta<-(mu[2]-mu[1])*2/3
  D<-diag(1,K-1)[1:(K-3),]+cbind(matrix(0,K-3,1),diag(-2,K-2)[1:(K-3),])+cbind(matrix(0,K-3,2),diag(1,K-3))
  P<-t(D)%*%D
  Ps<-P+diag(c(1,1,rep(0,K-3)))
  # Gibbs采样
  nu<-2
  A<-10^5
  z<-matrix(1,ncol=n,nrow=iter)
  a<-rep(1,iter)
  tau<-rep(1,iter)
  beta<-matrix(1,ncol=K-1,nrow=iter)
  for (i in 2:iter) {
    c<-exp(c(0,beta[i-1,]))/sum(exp(c(0,beta[i-1,])))
    for (p in 1:n) {
      H<-c*exp(-(x[p]-mu)^2/2/delta^2)/sum(c*exp(-(x[p]-mu)^2/2/delta^2))
      z[i,p]<-sample(1:K,1,replace=T,prob=H)
    }
    a[i]<-1/rgamma(1,shape=(nu+1)/2,rate=1/A^2+nu/tau[i-1])
    tau[i]<-1/rgamma(1,shape=(K+nu-1)/2,rate=nu/a[i]+beta[i-1,]%*%Ps%*%beta[i-1,]/2)
    beta[i,]<-hmc(tau[i],z[i,],beta[i-1,],K,Ps,n)
  }
  b<-burnin+1
  beta.hat<-colMeans(beta[b:iter,])
  # 模型代入
  c.hat<-exp(c(0,beta.hat))/sum(exp(c(0,beta.hat)))
  x0<-seq(a0,b0,length.out=n)
  f.hat<-rep(0,n)
  for (i in 1:n) {
    f.hat[i]<-sum(c.hat*exp(-(x0[i]-mu)^2/2/delta^2)/sqrt(2*pi*delta^2))
  }
  # 返回值
  return(list(t=x0,pgm=f.hat))
}
UU<-function(b,nn,c,tau0,Ps){
  return(-sum(nn*log(c))+t(b)%*%Ps%*%b/2/tau0)
}
KK<-function(u){
  return(t(u)%*%u/2)
}
hmc<-function(tau0,z0,beta0,K,Ps,n,v=0.018,L=10){
  ss<-K-1
  M<-diag(1,ss)
  u0<-matrix(rmvnorm(1,mu=rep(0,ss),sigma = M)) 
  #beta0<-matrix(1,K-1)
  nn<-rep(0,K)
  for (j in 1:K) {
    nn[j]<-sum(z0==j)
  }
  u1<-u0
  beta1<-beta0
  c<-exp(c(0,beta1))/sum(exp(c(0,beta1)))
  duc<-nn[-1]-n*c[-1]
  for (l in 1:L) {
    u1<-u1-v/2*(-duc+Ps%*%beta1/tau0)
    beta1<-beta1+v*u1
    u1<-u1-v/2*(-duc+Ps%*%beta1/tau0)
  }   
  p<-min(1,exp(-UU(beta1,nn,c,tau0,Ps)+UU(beta0,nn,c,tau0,Ps)-KK(u1)+KK(u0)))
  #if(rbinom(1,1,p)==1){
  if(p>0){
    return(beta1)
  }else{
    return(beta0)
  }
}