fashr/check_package.Rmd at main · stephenslab/fashr · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
---
title: 'Simulation: Pollution (Package)'
output: html_document
date: "2024-12-28"
---

```{r}
library(BayesGP)
library(TMB)
library(Matrix)
library(splines)
library(parallel)
library(ggplot2)
library(reshape2)
library(mixsqp)
library(tidyverse)
library(fashr)
function_dir <- paste0("/Users/ziangzhang/Desktop/FASH/FASHresultsummary", "/code/function")
source(paste0(function_dir, "/functions_fitting_Gaussian.R"))
source(paste0(function_dir, "/functions_simulation.R"))
```

```{r}
num_cores <- 1
# Generate indices of groups
group_indices <- rep(1:3, each = 16)

# Generate datasets in parallel
datasets <- mclapply(1:48, function(i) {
  set.seed(i)
  if (i <= 16) {
    n_basis <- 5
    sd_fun <- 1
  } else if (i <= 32) {
    n_basis <- 10
    sd_fun <- 1
  } else {
    n_basis <- 20
    sd_fun <- 1
  }
  simulate_process(n = 100, n_basis = n_basis, sd_fun = sd_fun, sd = 0.1)
}, mc.cores = num_cores)
```


```{r}
set.seed(123)
p_vec <- 2
log_prec <- unique(sort(c(Inf, seq(-1,1, by = 0.1), seq(-5,-1, by = 0.5), seq(1,5, by = 0.5), seq(-10,-5, by = 1), seq(5,10, by = 1)), decreasing = T))
psd_iwp_vec <- 1/exp(.5*log_prec)
```


```{r}
# using fashr
fashr_result <- fashr(data_list = datasets, Y = "y", smooth_var = "x",
                      S = 0.1, likelihood = "gaussian", order = 2,
                      grid = psd_iwp_vec, betaprec = 0.001, num_basis = 50,
                      verbose = TRUE)
plot(fashr_result)
```


```{r}
i <- 1
predi <- predict(fashr_result, index = i)
plot(datasets[[i]]$x, datasets[[i]]$y, type = 'p',
       main = paste("Dataset", i), xlab = "x", ylab = "y",
       cex = 0.1,
       cex.main = 0.8, cex.lab = 0.7, cex.axis = 0.7)
lines(datasets[[i]]$x, datasets[[i]]$truef, col = "red")
lines(predi$x, predi$mean, col = "blue")
lines(predi$x, predi$lower, col = "blue", lty = 2)
lines(predi$x, predi$upper, col = "blue", lty = 2)
```


```{r}
i <- 1
predi_samps <- predict(fashr_result, index = i, only.samples = TRUE, M = 5)
plot(datasets[[i]]$x, datasets[[i]]$y, type = 'p',
       main = paste("Dataset", i), xlab = "x", ylab = "y",
       cex = 0.1,
       cex.main = 0.8, cex.lab = 0.7, cex.axis = 0.7)
lines(datasets[[i]]$x, datasets[[i]]$truef, col = "red")
matlines(datasets[[i]]$x, predi_samps, col = "blue", lty = 1, lwd = 0.3)
```


```{r}
i <- 1
refined_x <- seq(from = 0, to = 5, length.out = 200)
predi_refined <- predict(fashr_result, index = i, smooth_var = refined_x)
plot(datasets[[i]]$x, datasets[[i]]$y, type = 'p',
       main = paste("Dataset", i), xlab = "x", ylab = "y",
       cex = 0.1,
       cex.main = 0.8, cex.lab = 0.7, cex.axis = 0.7)
lines(datasets[[i]]$x, datasets[[i]]$truef, col = "red")
lines(predi_refined$x, predi_refined$mean, col = "blue")
lines(predi_refined$x, predi_refined$lower, col = "blue", lty = 2)
lines(predi_refined$x, predi_refined$upper, col = "blue", lty = 2)
```