20250819 meeting #2

 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

---
title: 'Re-make the simulation study in Morris et al. (2017)'
output: html_notebook
date: "2025-04-21"
---

# Reference

Samuel A. Morris, Brian J. Reich, Emeric Thibaud, Daniel Cooley, A Space-Time Skew-t Model for Threshold Exceedances, Biometrics, Volume 73, Issue 3, September 2017, Pages 749–758, https://doi.org/10.1111/biom.12644

# Library and source
# libraries
library(fields)
library(SpatialTools)

# necessary functions
source("mcmc.R")
source("auxfunctions.R")

set.seed(100)
ns <- 2331
nt <- 260
s <- cbind(runif(n = ns, min = 0, max = 10), runif(n = ns, min = 0, max = 10))
x <- array(1, c(ns, nt, 3))
for (t in 1:nt) {
  x[, t, 2] <- s[, 1]
  x[, t, 3] <- s[, 2]
}
plot(s)


# generated data
set.seed(123)

lambda    <- c(0, 0,  0,  0,  3,  3,  3,  3,   0,   0,   0,   0,   3,   3,   3,   3,   0,   0,   0,   0,   3,   3,   3,   3)
tau.alpha <- c(0, 0,  0,  0,  6,  6,  6,  6,   0,   0,   0,   0,   6,   6,   6,   6,   0,   0,   0,   0,   6,   6,   6,   6)
tau.beta  <- c(0, 0,  0,  0, 16, 16, 16, 16,   0,   0,   0,   0,  16,  16,  16,  16,   0,   0,   0,   0,  16,  16,  16,  16)
nknots    <- c(1, 5, 10, 20,  1,  5, 10, 20,   1,   5,  10,  20,   1,   5,  10,  20,   1,   5,  10,  20,   1,   5,  10,  20)
phi.z     <- c(0, 0,  0,  0,  0,  0,  0,  0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9)
phi.w     <- c(0, 0,  0,  0,  0,  0,  0,  0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9)
phi.tau   <- c(0, 0,  0,  0,  0,  0,  0,  0, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5)

combined = array(NA, c(ns, nt, 24))

for (i in 1:24) {
  data_ <- rpotspatTS(nt = nt,
                      x = x, 
                      s = s, 
                      beta = c(10, 0, 0), 
                      gamma = 0.9, 
                      nu = 0.5, 
                      rho = 1, 
                      phi.z = phi.z[i],      # time series parameters
                      phi.w = phi.w[i],      # time series parameters
                      phi.tau = phi.tau[i],  # time series parameters
                      lambda = lambda[i], 
                      tau.alpha = tau.alpha[i], 
                      tau.beta = tau.beta[i], 
                      nknots = nknots[i], 
                      dist = "gaussian"
                      )$y
  
  # combined
  combined[, , i] = data_
  
  # hist
  hist(data_, main = paste0('hist. of ', i))
  
  # remove
  rm(data_)
}

# combine
#library(abind)
#combined <- abind(data_1, data_2, data_3, along = 3)

# save
library(dplyr)
library(reticulate)
np = import("numpy")
combined = combined %>% r_to_py()
combined = combined$transpose(2L, 1L, 0L)
save_dir = 'datasets'
if (!dir.exists(save_dir)) {
  dir.create(save_dir, recursive = TRUE)
}
np$save(paste0(save_dir, '/real.npy'), combined)
np$save(paste0(save_dir, '/all_locations.npy'), s)