doc/intro.R

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

## ----eval=FALSE---------------------------------------------------------------
#  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))
#  }

## ----eval=FALSE---------------------------------------------------------------
#  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)
#    }
#  }