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

Overview

StatComp21080 is a simple R package developed to estimate the density by using Lindsey's method and penalized Gaussian mixtures for the 'Statistical Computing' course. Two functions are considered, namely, LM (density estimation using Lindsey's method) and PGM (Density estimations using Lindsey's method and penalized Gaussian mixtures). For each function, Gibbs sampling is uesd.

function LM

The source R code for LM is as follows:

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

function PGM

The source R code for PGM is as follows:

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

SunJinglan/StatComp21080 documentation built on Dec. 24, 2021, 1:24 a.m.