docs/HMTree_demo.R

In shimlab/HMTree: Hidden Markov tree prior for analysis of multiple functional data

#Same as .Rmd (at least should be similar)

#install.packages("devtools")
#install.packages("Rcpp")
#install_github("shimlab/HMTree")
library("devtools")
library("Rcpp")
library("HMTree")

#Simulate Data
curve.length = 1024
model.mode = 'Poisson'
num.samples=1
set.seed(666)
res.pois = simu.curves(curve.length=curve.length, model.mode=model.mode, num.samples=num.samples)

curve.length = 1024
model.mode = 'Normal'
num.samples=1
normal.sigma=0.5
set.seed(666)
res.normal = simu.curves(curve.length=curve.length, model.mode=model.mode, num.samples=num.samples, normal.sigma=normal.sigma)

#After simulating data, you can get denoised objects by:

denoised.normal<-HMTree.denoise_1D(res.normal)
denoised.pois<-HMTree.denoise_1D(res.pois)


local.path<-("/Users/haosicheng/Desktop/") #Only need to change "local.path"
path<-paste0(local.path,"HMTree/docs/output")
#Normal:
#wavelet coef
path_normal_wavelet_coef<-paste0(path,"/hmt_1024_normal.wavelets_coef.txt")
normal_wavelet_coef <- read.table(path_normal_wavelet_coef,col.names = 1, quote="\"", comment.char="")
#denoised wavelet coef
path_normal_denoised_wavelet_coef<-paste0(path,"/hmt_1024_normal_denoised.wavelets_coef.txt")
normal_denoised_wavelet_coef <- read.table(path_normal_denoised_wavelet_coef,col.names = 1, quote="\"", comment.char="")
#denoised data
path_normal_denoised<-paste0(path,"/hmt_1024_normal_denoised.txt")
normal_denoised <- read.table(path_normal_denoised,col.names = 1, quote="\"", comment.char="")
#Poisson:
#wavelet coef
path_pois_wavelet_coef<-paste0(path,"/pois_1024.wavelets_coef.txt")
pois_wavelet_coef<-read.table(path_pois_wavelet_coef,col.names = 1, quote="\"", comment.char="")
#denoised wavelet coef
path_pois_denoised_wavelet_coef<-paste0(path,"/pois_1024_denoised.wavelets_coef.txt")
pois_denoised_wavelet_coef <- read.table(path_pois_denoised_wavelet_coef,col.names = 1, quote="\"", comment.char="")
#denoised data
path_pois_denoised<-paste0(path,"/pois_1024_denoised.txt")
pois_denoised <- read.table(path_pois_denoised,col.names = 1, quote="\"", comment.char="")


#Compare:

#Normal
normal.output<-denoised.normal
diff.normal.wavelet.coef<-(normal.output$wavelet_coef-normal_wavelet_coef)/normal_wavelet_coef
summary(diff.normal.wavelet.coef)

diff.normal.denoised.wavelet.coef<-(normal.output$denoised_wavelet_coef-normal_denoised_wavelet_coef)/normal_denoised_wavelet_coef
summary(diff.normal.wavelet.coef)

diff.normal.denoised<-(normal.output$denoised-normal_denoised)/normal_denoised
summary(diff.normal.denoised)
#Poisson
pois.output<-denoised.pois
diff.pois.wavelet.coef<-(pois.output$wavelet_coef-pois_wavelet_coef)/pois_wavelet_coef
summary(diff.pois.wavelet.coef)

diff.pois.denoised.wavelet.coef<-(pois.output$denoised_wavelet_coef-pois_denoised_wavelet_coef)/pois_denoised_wavelet_coef
summary(diff.pois.wavelet.coef)

diff.pois.denoised<-(pois.output$denoised-pois_denoised)/pois_denoised
summary(diff.pois.denoised)