output_4_ls.Rout

R version 4.2.3 (2023-03-15) -- "Shortstop Beagle"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: aarch64-apple-darwin20 (64-bit)

R é um software livre e vem sem GARANTIA ALGUMA.
Você pode redistribuí-lo sob certas circunstâncias.
Digite 'license()' ou 'licence()' para detalhes de distribuição.

R é um projeto colaborativo com muitos contribuidores.
Digite 'contributors()' para obter mais informações e
'citation()' para saber como citar o R ou pacotes do R em publicações.

Digite 'demo()' para demonstrações, 'help()' para o sistema on-line de ajuda,
ou 'help.start()' para abrir o sistema de ajuda em HTML no seu navegador.
Digite 'q()' para sair do R.

> ##### Oliver's MAC ####
> reticulate::use_python("/Users/oli/miniforge3/envs/r-tensorflow/bin/python3.8", required = TRUE)
> library(reticulate)
> reticulate::py_config()
python:         /Users/oli/miniforge3/envs/r-tensorflow/bin/python3.8
libpython:      /Users/oli/miniforge3/envs/r-tensorflow/lib/libpython3.8.dylib
pythonhome:     /Users/oli/miniforge3/envs/r-tensorflow:/Users/oli/miniforge3/envs/r-tensorflow
version:        3.8.18 (default, Sep 11 2023, 08:17:16)  [Clang 14.0.6 ]
numpy:          /Users/oli/miniforge3/envs/r-tensorflow/lib/python3.8/site-packages/numpy
numpy_version:  1.24.3

NOTE: Python version was forced by use_python() function
> 
> # Get command-line arguments - if called via sh
> args <- commandArgs(trailingOnly = TRUE)
> if (length(args) == 0) {  # if not called via sh
+   args <- c(5, 'cs') 
+ }
> F32 <- as.numeric(args[1])
> M32 <- args[2]
> print(paste("F32:", F32, "M32:", M32))
[1] "F32: 5 M32: ls"
> 
> 
> #### A mixture of discrete and continuous variables ####
> library(tensorflow)
> library(keras)
> library(mlt)
Carregando pacotes exigidos: basefun
Carregando pacotes exigidos: variables
> library(tram)
> library(MASS)
> library(tensorflow)
> library(keras)
> library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.1     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::desc()   masks variables::desc()
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
✖ dplyr::select() masks MASS::select()
✖ ggplot2::unit() masks variables::unit()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
> source('summerof24/utils_tf.R')

Attaching package: ‘grid’

The following object is masked from ‘package:variables’:

    unit

> 
> #### For TFP
> library(tfprobability)
> source('summerof24/utils_tfp.R')
> 
> ##### Flavor of experiment ######
> 
> #### Saving the current version of the script into runtime
> DIR = 'summerof24/runs/triangle_structured_continous/run_nodes25'
> if (!dir.exists(DIR)) {
+   dir.create(DIR, recursive = TRUE)
+ }
> # Copy this file to the directory DIR
> file.copy('summerof24/triangle_structured_continous.R', file.path(DIR, 'triangle_structured_continous.R'), overwrite=TRUE)
[1] TRUE
> 
> num_epochs <- 500
> len_theta = 20 # Number of coefficients of the Bernstein polynomials
> hidden_features_I = c(2,25,25,2)    #hidden_features_CS=hidden_features_I = c(2,25,25,2)
> hidden_features_CS = c(2,25,25,2)
> 
> SEED = -1 #If seed > 0 then the seed is set
> 
> if (F32 == 1){
+   FUN_NAME = 'DPGLinear'
+   f <- function(x) -0.3 * x
+ } else if (F32 == 2){
+   f = function(x) 2 * x**3 + x
+   FUN_NAME = 'DPG2x3+x'
+ } else if (F32 == 3){
+   f = function(x) 0.5*exp(x)
+   FUN_NAME = 'DPG0.5exp'
+ } else if (F32 == 4){
+   f = function(x) 0.75*atan(5*(x+0.12)) 
+   FUN_NAME = 'DPGatan'
+ } else if (F32 == 5){
+   f = function(x) 2*sin(3*x)+x 
+   FUN_NAME = 'DPGSin'
+ } else {
+   stop("Unknown Function F32")
+ }
> 
> # xs = seq(-1,1,0.1)
> # f = function(x) 2 * x**3 + 0.1*(x-.5)**5
> if (FALSE){
+   f = function(x) 0.75*atan(5*(x+0.12)) 
+   plot(xs, f(xs))
+   s=train$df_R$x2
+   hist(s, freq=FALSE, 100)
+   hist(f(s), freq=FALSE, 100)
+ }
> 
> if (M32 == 'ls') {
+   MA =  matrix(c(
+     0, 'ls', 'ls', 
+     0,    0, 'ls', 
+     0,    0,   0), nrow = 3, ncol = 3, byrow = TRUE)
+   MODEL_NAME = 'ModelLS'
+ } else{
+   MA =  matrix(c(
+     0, 'ls', 'ls', 
+     0,    0, 'cs', 
+     0,    0,   0), nrow = 3, ncol = 3, byrow = TRUE)
+   MODEL_NAME = 'ModelCS'
+ }
> 
> 
> # fn = 'triangle_mixed_DGPLinear_ModelLinear.h5'
> # fn = 'triangle_mixed_DGPSin_ModelCS.h5'
> 
> if (SEED < 0){
+   fn = file.path(DIR, paste0('triangle_mixed_', FUN_NAME, '_', MODEL_NAME))
+ } else{
+   fn = file.path(DIR, paste0('triangle_mixed_', FUN_NAME, '_', MODEL_NAME, '_SEED', SEED))
+ }
> print(paste0("Starting experiment ", fn))
[1] "Starting experiment summerof24/runs/triangle_structured_continous/run_nodes25/triangle_mixed_DPGSin_ModelLS"
>    
> xs = seq(-1,1,0.1)
> 
> plot(xs, -f(xs), sub=fn, xlab='x2', ylab='f(x2)', main='DGP influence of x2 on x3', cex.sub=0.4)
> ##### DGP ########
> dgp <- function(n_obs, doX=c(NA, NA, NA), seed=-1) {
+     if (seed > 0) {
+       set.seed(seed)
+       print(paste0("Setting Seed:", seed))
+     }
+     #n_obs = 1e5 n_obs = 10
+     #Sample X_1 from GMM with 2 components
+     if (is.na(doX[1])){
+       X_1_A = rnorm(n_obs, 0.25, 0.1)
+       X_1_B = rnorm(n_obs, 0.73, 0.05)
+       X_1 = ifelse(sample(1:2, replace = TRUE, size = n_obs) == 1, X_1_A, X_1_B)
+     } else{
+       X_1 = rep(doX[1], n_obs)
+     }
+     #hist(X_1)
+     
+     # Sampling according to colr
+     if (is.na(doX[2])){
+       U2 = runif(n_obs)
+       
+       x_2_dash = qlogis(U2)
+       #x_2_dash = h_0(x_2) + beta * X_1
+       X_2 = 1/0.42 * (x_2_dash - 2 * X_1)
+       X_2 = 1/5. * (x_2_dash - 0.4 * X_1) # 0.39450
+       X_2 = 1/5. * (x_2_dash - 1.2 * X_1) 
+       #h2 = x_2_dash = 5 * x_2 + 2 * X_1
+       X_2 = 1/5. * (x_2_dash - 2 * X_1)  # 
+       
+       
+     } else{
+       X_2 = rep(doX[2], n_obs)
+     }
+     
+     #hist(X_2)
+     #ds = seq(-5,5,0.1)
+     #plot(ds, dlogis(ds))
+     if (is.na(doX[3])){
+       U3 = runif(n_obs)
+       x_3_dash = qlogis(U3)
+       #h(x3|x1,x2) = 0.63*x3 - 0.2*x1 - f(x2)
+       #x_3_dash = h_0_3(x_3) + gamma_1 * X_1 + gamma_2 * X_2
+       #x_3_dash = 0.63 * x_3 -0.2 * X_1 + 1.3 * X_2
+       #x_3_dash = h(x3|x1,x2) = 0.63*x3 - 0.2*x1 - f(x2)
+       X_3 = (x_3_dash + 0.2 * X_1 + f(X_2))/0.63
+     } else{
+       X_3 = rep(doX[3], n_obs)
+     }
+     
+    
+     #hist(X_3)
+     A <- matrix(c(0, 1, 1, 0,0,1,0,0,0), nrow = 3, ncol = 3, byrow = TRUE)
+     dat.orig =  data.frame(x1 = X_1, x2 = X_2, x3 = X_3)
+     dat.tf = tf$constant(as.matrix(dat.orig), dtype = 'float32')
+     
+     q1 = quantile(dat.orig[,1], probs = c(0.05, 0.95)) 
+     q2 = quantile(dat.orig[,2], probs = c(0.05, 0.95))
+     q3 = quantile(dat.orig[,3], probs = c(0.05, 0.95))
+     
+     
+     return(list(
+       df_orig=dat.tf, 
+       df_R = dat.orig,
+       #min =  tf$reduce_min(dat.tf, axis=0L),
+       #max =  tf$reduce_max(dat.tf, axis=0L),
+       min = tf$constant(c(q1[1], q2[1], q3[1]), dtype = 'float32'),
+       max = tf$constant(c(q1[2], q2[2], q3[2]), dtype = 'float32'),
+       type = c('c', 'c', 'c'),
+       A=A))
+ } 
> 
> train = dgp(40000, seed=ifelse(SEED > 0, SEED, -1))
> test  = dgp(40000, seed=ifelse(SEED > 0, SEED + 1, -1))
> (global_min = train$min)
tf.Tensor([ 0.12299408 -0.8056552  -6.966295  ], shape=(3), dtype=float32)
> (global_max = train$max)
tf.Tensor([0.7941454  0.41015792 5.2153287 ], shape=(3), dtype=float32)
> data_type = train$type
> 
> 
> 
> #### Fitting Tram ######
> df = data.frame(train$df_orig$numpy())
> fit.orig = Colr(X2~X1, order=len_theta ,df)
> summary(fit.orig)

  Continuous Outcome Logistic Regression 

Call:
Colr(formula = X2 ~ X1, data = df, order = len_theta)

Coefficients:
   Estimate Std. Error z value Pr(>|z|)    
X1  2.03408    0.03543   57.41   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log-Likelihood:
 -15358.07 (df = 22)
Likelihood-ratio Test: Chisq = 3376.035 on 1 degrees of freedom; p = < 2.2e-16

> confint(fit.orig) #Original
      2.5 %  97.5 %
X1 1.964636 2.10353
> #dd = predict(fit.orig, newdata = data.frame(X1 = 0.5), type = 'density')
> #x2s = as.numeric(rownames(dd))
> #plot(x2s, dd, type = 'l', col='red')
> 
> #?predict.tram
> summary(fit.orig)

  Continuous Outcome Logistic Regression 

Call:
Colr(formula = X2 ~ X1, data = df, order = len_theta)

Coefficients:
   Estimate Std. Error z value Pr(>|z|)    
X1  2.03408    0.03543   57.41   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log-Likelihood:
 -15358.07 (df = 22)
Likelihood-ratio Test: Chisq = 3376.035 on 1 degrees of freedom; p = < 2.2e-16

> confint(fit.orig) #Original 
      2.5 %  97.5 %
X1 1.964636 2.10353
> 
> # Fitting Tram
> fit.orig = Colr(x3 ~ x1 + x2 ,order=len_theta ,train$df_R)
> summary(fit.orig)

  Continuous Outcome Logistic Regression 

Call:
Colr(formula = x3 ~ x1 + x2, data = train$df_R, order = len_theta)

Coefficients:
   Estimate Std. Error z value Pr(>|z|)    
x1 -0.01420    0.03546    -0.4    0.689    
x2 -3.17932    0.02992  -106.3   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log-Likelihood:
 -102543.9 (df = 23)
Likelihood-ratio Test: Chisq = 13834.98 on 2 degrees of freedom; p = < 2.2e-16

> confint(fit.orig) #Original 
         2.5 %      97.5 %
x1 -0.08369958  0.05530882
x2 -3.23796316 -3.12067455
> 
> len_theta_max = len_theta
> for (i in 1:nrow(MA)){ #Maximum number of coefficients (BS and Levels - 1 for the ordinal)
+   if (train$type[i] == 'o'){
+     len_theta_max = max(len_theta_max, nlevels(train$df_R[,i]) - 1)
+   }
+ }
> param_model = create_param_model(MA, hidden_features_I = hidden_features_I, 
+                                  len_theta = len_theta_max, 
+                                  hidden_features_CS = hidden_features_CS)
> 
> x = tf$ones(shape = c(2L, 3L))
> param_model(1*x)
tf.Tensor(
[[[ 0.00000000e+00  0.00000000e+00 -4.09818925e-02 -4.82878201e-02
   -2.03795340e-02  2.04079691e-02 -1.45959547e-02  5.15200682e-02
   -1.76881216e-02 -2.70995703e-02  1.05715748e-02  1.47924917e-02
   -1.04625747e-02  9.54719558e-02  1.59317359e-01  3.52589898e-02
   -1.48528991e-02 -1.06845830e-04 -1.34109741e-03  6.23627491e-02
    4.51825261e-02  2.21384540e-02]
  [ 0.00000000e+00 -1.07005395e-01 -1.51259517e-02  5.51365204e-02
   -1.71946350e-03 -2.71369126e-02 -4.56792787e-02 -4.68725991e-03
    4.70794598e-03  1.00934086e-02  1.47555336e-01 -2.59636398e-02
   -7.71734267e-02 -1.34865120e-02 -1.01983741e-01 -1.14899715e-02
    1.98853128e-02  8.39964822e-02 -6.65639266e-02 -4.44038212e-02
    9.40392017e-02  1.76062267e-02]
  [ 0.00000000e+00  8.46536271e-03  2.52183876e-03 -2.26724963e-03
   -5.59973791e-02  8.33005272e-03  7.44136199e-02 -2.76201703e-02
    3.71399596e-02  5.45193069e-02 -3.85128148e-03  5.09912148e-02
   -3.49550955e-02 -6.43406436e-02  6.40737936e-02  7.51333265e-03
   -1.36419032e-02  1.46925328e-02  4.98599792e-03  3.40193808e-02
   -2.65461709e-02 -2.85198838e-02]]

 [[ 0.00000000e+00  0.00000000e+00 -4.09818925e-02 -4.82878201e-02
   -2.03795340e-02  2.04079691e-02 -1.45959547e-02  5.15200682e-02
   -1.76881216e-02 -2.70995703e-02  1.05715748e-02  1.47924917e-02
   -1.04625747e-02  9.54719558e-02  1.59317359e-01  3.52589898e-02
   -1.48528991e-02 -1.06845830e-04 -1.34109741e-03  6.23627491e-02
    4.51825261e-02  2.21384540e-02]
  [ 0.00000000e+00 -1.07005395e-01 -1.51259517e-02  5.51365204e-02
   -1.71946350e-03 -2.71369126e-02 -4.56792787e-02 -4.68725991e-03
    4.70794598e-03  1.00934086e-02  1.47555336e-01 -2.59636398e-02
   -7.71734267e-02 -1.34865120e-02 -1.01983741e-01 -1.14899715e-02
    1.98853128e-02  8.39964822e-02 -6.65639266e-02 -4.44038212e-02
    9.40392017e-02  1.76062267e-02]
  [ 0.00000000e+00  8.46536271e-03  2.52183876e-03 -2.26724963e-03
   -5.59973791e-02  8.33005272e-03  7.44136199e-02 -2.76201703e-02
    3.71399596e-02  5.45193069e-02 -3.85128148e-03  5.09912148e-02
   -3.49550955e-02 -6.43406436e-02  6.40737936e-02  7.51333265e-03
   -1.36419032e-02  1.46925328e-02  4.98599792e-03  3.40193808e-02
   -2.65461709e-02 -2.85198838e-02]]], shape=(2, 3, 22), dtype=float32)
> MA
     [,1] [,2] [,3]
[1,] "0"  "ls" "ls"
[2,] "0"  "0"  "ls"
[3,] "0"  "0"  "0" 
> h_params = param_model(train$df_orig)
> # Check the derivatives of h w.r.t. x
> x <- tf$ones(shape = c(2L, 3L)) #B,P
> with(tf$GradientTape(persistent = TRUE) %as% tape, {
+   tape$watch(x)
+   y <- param_model(x)
+ })
> # parameter (output) has shape B, P, k (num param)
> # derivation of param wrt to input x
> # input x has shape B, P
> # derivation d has shape B,P,k, B,P
> d <- tape$jacobian(y, x)
> d[1,,,2,] # only contains zero since independence of batches
tf.Tensor(
[[[0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]]

 [[0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]]

 [[0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]
  [0. 0. 0.]]], shape=(3, 22, 3), dtype=float32)
> 
> 
> # loss before training
> struct_dag_loss(t_i=train$df_orig, h_params=h_params)
tf.Tensor(2.5264359, shape=(), dtype=float32)
Warning message:
In force(if_any_TRUE) :
  Indexing tensors are passed as-is to python, no index offsetting or R to python translation is performed. Selected options for one_based and inclusive_stop are ignored and treated as FALSE. To silence this warning, set options(tensorflow.extract.warn_tensors_passed_asis = FALSE)
> 
> if (FALSE){
+   with(tf$GradientTape(persistent = TRUE) %as% tape, {
+     h_params = param_model(train$df_orig)
+     loss = struct_dag_loss(train$df_orig, h_params)
+   })
+   gradients = tape$gradient(loss, param_model$trainable_variables)
+   gradients
+ }
> 
> param_model = create_param_model(MA, hidden_features_I=hidden_features_I, len_theta=len_theta, hidden_features_CS=hidden_features_CS)
> 
> 
> # ######### DEBUG TRAINING FROM HAND #######
> # # Define the optimizer
> # optimizer <- tf$optimizers$Adam(lr=0.01)
> # # Define the number of epochs for training
> # num_epochs <- 10
> # for (epoch in 1:num_epochs) {
> #   with(tf$GradientTape(persistent = TRUE) %as% tape, {
> #     # Compute the model's prediction - forward pass
> #     h_params <- param_model(train$df_scaled)
> #     loss <- struct_dag_loss(train$df_scaled, h_params)
> #   })
> #   # Compute gradients
> #   gradients <- tape$gradient(loss, param_model$trainable_variables)
> #   # Apply gradients to update the model parameters
> #   optimizer$apply_gradients(purrr::transpose(list(gradients, param_model$trainable_variables)))
> #   # Print the loss every epoch or more frequently if desired
> #   print(paste("Epoch", epoch, ", Loss:", loss$numpy()))
> # }
> 
> 
> optimizer = optimizer_adam()
WARNING:absl:At this time, the v2.11+ optimizer `tf.keras.optimizers.Adam` runs slowly on M1/M2 Macs, please use the legacy Keras optimizer instead, located at `tf.keras.optimizers.legacy.Adam`.
> param_model$compile(optimizer, loss=struct_dag_loss)
WARNING:absl:There is a known slowdown when using v2.11+ Keras optimizers on M1/M2 Macs. Falling back to the legacy Keras optimizer, i.e., `tf.keras.optimizers.legacy.Adam`.
> param_model$evaluate(x = train$df_orig, y=train$df_orig, batch_size = 7L)

   1/5715 [..............................] - ETA: 36:45 - loss: 2.5282
  88/5715 [..............................] - ETA: 3s - loss: 2.5335   
 216/5715 [>.............................] - ETA: 2s - loss: 2.5346
 345/5715 [>.............................] - ETA: 2s - loss: 2.5346
 475/5715 [=>............................] - ETA: 2s - loss: 2.5347
 605/5715 [==>...........................] - ETA: 2s - loss: 2.5348
 737/5715 [==>...........................] - ETA: 2s - loss: 2.5348
 868/5715 [===>..........................] - ETA: 1s - loss: 2.5348
 998/5715 [====>.........................] - ETA: 1s - loss: 2.5348
1126/5715 [====>.........................] - ETA: 1s - loss: 2.5347
1256/5715 [=====>........................] - ETA: 1s - loss: 2.5347
1386/5715 [======>.......................] - ETA: 1s - loss: 2.5348
1513/5715 [======>.......................] - ETA: 1s - loss: 2.5347
1644/5715 [=======>......................] - ETA: 1s - loss: 2.5348
1775/5715 [========>.....................] - ETA: 1s - loss: 2.5348
1904/5715 [========>.....................] - ETA: 1s - loss: 2.5348
2034/5715 [=========>....................] - ETA: 1s - loss: 2.5348
2164/5715 [==========>...................] - ETA: 1s - loss: 2.5348
2300/5715 [===========>..................] - ETA: 1s - loss: 2.5348
2431/5715 [===========>..................] - ETA: 1s - loss: 2.5348
2561/5715 [============>.................] - ETA: 1s - loss: 2.5348
2692/5715 [=============>................] - ETA: 1s - loss: 2.5348
2822/5715 [=============>................] - ETA: 1s - loss: 2.5348
2954/5715 [==============>...............] - ETA: 1s - loss: 2.5348
3085/5715 [===============>..............] - ETA: 1s - loss: 2.5348
3216/5715 [===============>..............] - ETA: 0s - loss: 2.5349
3346/5715 [================>.............] - ETA: 0s - loss: 2.5348
3475/5715 [=================>............] - ETA: 0s - loss: 2.5348
3606/5715 [=================>............] - ETA: 0s - loss: 2.5349
3737/5715 [==================>...........] - ETA: 0s - loss: 2.5349
3869/5715 [===================>..........] - ETA: 0s - loss: 2.5349
3999/5715 [===================>..........] - ETA: 0s - loss: 2.5349
4129/5715 [====================>.........] - ETA: 0s - loss: 2.5349
4258/5715 [=====================>........] - ETA: 0s - loss: 2.5349
4389/5715 [======================>.......] - ETA: 0s - loss: 2.5349
4520/5715 [======================>.......] - ETA: 0s - loss: 2.5349
4651/5715 [=======================>......] - ETA: 0s - loss: 2.5349
4782/5715 [========================>.....] - ETA: 0s - loss: 2.5349
4915/5715 [========================>.....] - ETA: 0s - loss: 2.5349
5045/5715 [=========================>....] - ETA: 0s - loss: 2.5348
5174/5715 [==========================>...] - ETA: 0s - loss: 2.5349
5304/5715 [==========================>...] - ETA: 0s - loss: 2.5348
5434/5715 [===========================>..] - ETA: 0s - loss: 2.5348
5566/5715 [============================>.] - ETA: 0s - loss: 2.5348
5700/5715 [============================>.] - ETA: 0s - loss: 2.5348
5715/5715 [==============================] - 3s 388us/step - loss: 2.5348
[1] 2.534813
> 
> 
> ##### Training or readin of weights if h5 available ####
> fnh5 = paste0(fn, '_E', num_epochs, '.h5')
> fnRdata = paste0(fn, '_E', num_epochs, '.RData')
> if (file.exists(fnh5)){
+   param_model$load_weights(fnh5)
+   load(fnRdata) #Loading of the workspace causes trouble e.g. param_model is zero
+   # Quick Fix since loading global_min causes problem (no tensors as RDS)
+   (global_min = train$min)
+   (global_max = train$max)
+ } else {
+   if (FALSE){ ### Full Training w/o diagnostics
+     hist = param_model$fit(x = train$df_orig, y=train$df_orig, epochs = 200L,verbose = TRUE)
+     param_model$save_weights(fn)
+     plot(hist$epoch, hist$history$loss)
+     plot(hist$epoch, hist$history$loss, ylim=c(1.07, 1.2))
+   } else { ### Training with diagnostics
+     ws <- data.frame(w12 = numeric())
+     train_loss <- numeric()
+     val_loss <- numeric()
+     
+     # Training loop
+     for (e in 1:num_epochs) {
+       print(paste("Epoch", e))
+       hist <- param_model$fit(x = train$df_orig, y = train$df_orig, 
+                               epochs = 1L, verbose = TRUE, 
+                               validation_data = list(test$df_orig,test$df_orig))
+       
+       # Append losses to history
+       train_loss <- c(train_loss, hist$history$loss)
+       val_loss <- c(val_loss, hist$history$val_loss)
+       
+       # Extract specific weights
+       w <- param_model$get_layer(name = "beta")$get_weights()[[1]]
+       
+       ws <- rbind(ws, data.frame(w12 = w[1, 2], w13 = w[1, 3], w23 = w[2, 3]))
+     }
+     # Save the model
+     param_model$save_weights(fnh5)
+     save(train_loss, val_loss, train_loss, f, MA, len_theta,
+          hidden_features_I,
+          hidden_features_CS,
+          ws,
+          #global_min, global_max,
+          file = fnRdata)
+   }
+ }
tf.Tensor([0.7941454  0.41015792 5.2153287 ], shape=(3), dtype=float32)
> 
> ####### FINISHED TRAINING #####
> #pdf(paste0('loss_',fn,'.pdf'))
> epochs = length(train_loss)
> plot(1:length(train_loss), train_loss, type='l', main='Training (black: train, green: valid)')
> lines(1:length(train_loss), val_loss, type = 'l', col = 'green')
> 
> # Last 50
> diff = max(epochs - 50,1)
> plot(diff:epochs, val_loss[diff:epochs], type = 'l', col = 'green', main='Last 50 epochs')
> lines(diff:epochs, train_loss[diff:epochs], type='l')
> 
> # plot(1:epochs, ws[,1], type='l', main='Coef', ylim=c(-0.5, 3))#, ylim=c(0, 6))
> # abline(h=2, col='green')
> # lines(1:epochs, ws[,2], type='l', ylim=c(0, 3))
> # abline(h=0.2, col='green')
> # lines(1:epochs, ws[,3], type='l', ylim=c(0, 3))
> # abline(h=-0.3, col='green')
> 
> 
> ggplot(ws, aes(x=1:nrow(ws))) + 
+   geom_line(aes(y=w12, color='x1 --> x2')) + 
+   geom_line(aes(y=w13, color='x1 --> x3')) + 
+   geom_line(aes(y=w23, color='x2 --> x3')) + 
+   geom_hline(aes(yintercept=2, color='x1 --> x2'), linetype=2) +
+   geom_hline(aes(yintercept=-0.2, color='x1 --> x3'), linetype=2) +
+   geom_hline(aes(yintercept=+0.3, color='x2 --> x3'), linetype=2) +
+   #scale_color_manual(values=c('x1 --> x2'='skyblue', 'x1 --> x3='red', 'x2 --> x3'='darkgreen')) +
+   labs(x='Epoch', y='Coefficients') +
+   theme_minimal() +
+   theme(legend.title = element_blank())  # Removes the legend title
> 
> if (FALSE){
+   p = ggplot(ws, aes(x=1:nrow(ws))) + 
+     geom_line(aes(y=w12, color="beta12")) + 
+     geom_line(aes(y=w13, color="beta13")) + 
+     #geom_line(aes(y=w23, color="beta23")) + 
+     geom_hline(aes(yintercept=2, color="beta12"), linetype=2) +
+     geom_hline(aes(yintercept=-0.2, color="beta13"), linetype=2) +
+     #geom_hline(aes(yintercept=+0.3, color="beta23"), linetype=2) +
+     scale_color_manual(
+       values=c('beta12'='skyblue', 'beta13'='red'),
+       labels=c(expression(beta[12]), expression(beta[13]))
+     ) +
+     labs(x='Epoch', y='Coefficients') +
+     theme_minimal() +
+     theme(
+       legend.title = element_blank(),   # Removes the legend title
+       legend.position = c(0.85, 0.25),  # Adjust this to position the legend inside the plot (lower-right)
+       legend.background = element_rect(fill="white", colour="black")  # Optional: white background with border
+     )
+   
+   file_name <- paste0(fn, "_coef_epoch.pdf")
+   # Save the plot
+   ggsave(file_name, plot = p, width = 8, height = 6)
+   file_path <- file.path("/Users/oli/Library/CloudStorage/Dropbox/Apps/Overleaf/tramdag/figures", basename(file_name))
+   ggsave(file_path, plot = p, width = 8/2, height = 6/2)
+ }
> 
> if (FALSE){
+ # Creating the figure for the paper 
+ # triangle_mixed_DPGLinear_ModelLS_coef_epoch 
+   p = ggplot(ws, aes(x=1:nrow(ws))) + 
+     geom_line(aes(y=w12, color="beta12")) + 
+     geom_line(aes(y=w13, color="beta13")) + 
+     geom_line(aes(y=w23, color="beta23")) + 
+     geom_hline(aes(yintercept=2, color="beta12"), linetype=2) +
+     geom_hline(aes(yintercept=-0.2, color="beta13"), linetype=2) +
+     geom_hline(aes(yintercept=+0.3, color="beta23"), linetype=2) +
+     scale_color_manual(
+       values=c('beta12'='skyblue', 'beta13'='red', 'beta23'='darkgreen'),
+       labels=c(expression(beta[12]), expression(beta[13]), expression(beta[23]))
+     ) +
+     labs(x='Epoch', y='Coefficients') +
+     theme_minimal() +
+     theme(
+       legend.title = element_blank(),   # Removes the legend title
+       legend.position = c(0.85, 0.25),  # Adjust this to position the legend inside the plot (lower-right)
+       legend.background = element_rect(fill="white", colour="black")  # Optional: white background with border
+     )
+   
+   file_name <- paste0(fn, "_coef_epoch.pdf")
+   # Save the plot
+   ggsave(file_name, plot = p, width = 8, height = 6)
+   file_path <- file.path("/Users/oli/Library/CloudStorage/Dropbox/Apps/Overleaf/tramdag/figures", basename(file_name))
+   ggsave(file_path, plot = p, width = 8/2, height = 6/2)
+ }
> 
> param_model$evaluate(x = train$df_orig, y=train$df_scaled) #Does not work, probably TF Eager vs Compiled

   1/1250 [..............................] - ETA: 1:29 - loss: 0.0000e+00
 170/1250 [===>..........................] - ETA: 0s - loss: 0.0000e+00  
 404/1250 [========>.....................] - ETA: 0s - loss: 0.0000e+00
 638/1250 [==============>...............] - ETA: 0s - loss: 0.0000e+00
 873/1250 [===================>..........] - ETA: 0s - loss: 0.0000e+00
1108/1250 [=========================>....] - ETA: 0s - loss: 0.0000e+00
1250/1250 [==============================] - 0s 225us/step - loss: 0.0000e+00
[1] 0
> # One more step to estimate NLL
> if (FALSE){
+   vals = NULL
+   for (i in 1:10){
+     test  = dgp(40000, i+10001)
+     hist = param_model$fit(x = train$df_orig, y = train$df_orig, 
+                     epochs = 1L, verbose = TRUE, 
+                     validation_data = list(test$df_orig,test$df_orig))
+     vals = append(vals, hist$history$val_loss)
+   }
+   t.test(vals)
+   M32
+   F32
+ }
> fn
[1] "summerof24/runs/triangle_structured_continous/run_nodes25/triangle_mixed_DPGSin_ModelLS"
> len_theta
[1] 20
> param_model$get_layer(name = "beta")$get_weights() * param_model$get_layer(name = "beta")$mask
tf.Tensor(
[[[ 0.          2.0012531   0.01728118]
  [ 0.          0.         -3.1162605 ]
  [ 0.          0.          0.        ]]], shape=(1, 3, 3), dtype=float32)
> 
> 
> #### Checking the transformation ####
> h_params = param_model(train$df_orig)
> r = check_baselinetrafo(h_params)
> Xs = r$Xs
> h_I = r$h_I
> 
> ##### X1
> fit.1 = Colr(X1~1,df, order=len_theta)
> plot(fit.1, which = 'baseline only', main='Black: COLR, Red: Our Model')
> lines(Xs[,1], h_I[,1], col='red', lty=2, lwd=3)
> rug(train$df_orig$numpy()[,1], col='blue')
> 
> 
> df = data.frame(train$df_orig$numpy())
> fit.21 = Colr(X2~X1,df, order=len_theta)
> temp = model.frame(fit.21)[1:2,-1, drop=FALSE] #WTF!
> plot(fit.21, which = 'baseline only', newdata = temp, lwd=2, col='blue', 
+      main='h_I(X2) Black: COLR, Red: Our Model', cex.main=0.8)
> lines(Xs[,2], h_I[,2], col='red', lty=2, lwd=5)
> rug(train$df_orig$numpy()[,2], col='blue')
> 
> fit.312 = Colr(X3 ~ X1 + X2,df, order=len_theta)
> temp = model.frame(fit.312)[1:2, -1, drop=FALSE] #WTF!
> 
> plot(fit.312, which = 'baseline only', newdata = temp, lwd=2, col='blue', 
+      main='h_I(X3) Colr and Our Model', cex.main=0.8)
> lines(Xs[,3], h_I[,3], col='red', lty=2, lwd=5)
> rug(train$df_orig$numpy()[,3], col='blue')
> 
> 
> if (FALSE){
+   # Check the derivatives of h w.r.t. x
+   x <- tf$ones(shape = c(10L, 3L)) #B,P
+   with(tf$GradientTape(persistent = TRUE) %as% tape, {
+     tape$watch(x)
+     y <- param_model(x)
+   })
+   d <- tape$jacobian(y, x)
+   for (k in 1:(2+len_theta)){ #k = 1
+     print(k) #B,P,k,B,P
+     B = 1
+     print(d[B,,k,B,]) #
+   }
+ }
> 
> ##### Checking observational distribution ####
> library(car)
Carregando pacotes exigidos: carData

Attaching package: ‘car’

The following object is masked from ‘package:dplyr’:

    recode

The following object is masked from ‘package:purrr’:

    some

> s = do_dag_struct(param_model, train$A, doX=c(NA, NA, NA), num_samples = 5000)
sample_from_target Fraction of extrapolated samples > 1 : %f 
0.196266666054726sample_from_target Fraction of extrapolated samples > 1 : %f 
0.213733330368996sample_from_target Fraction of extrapolated samples > 1 : %f 
0.0480666682124138> par(mfrow=c(1,3))
> for (i in 1:3){
+   d = s[,i]$numpy()
+   hist(train$df_orig$numpy()[,i], freq=FALSE, 100,main=paste0("X",i, " red: ours, black: data"), xlab='samples')
+   lines(density(train$df_orig$numpy()[,i]), col='blue')
+   #hist(train$df_orig$numpy()[,i], freq=FALSE, 100,main=paste0("X_",i))
+   lines(density(s[,i]$numpy()), col='red')
+   #qqplot(train$df_orig$numpy()[,i], s[,i]$numpy())
+   #abline(0,1)
+ }
> par(mfrow=c(1,1))
> 
> ######### Simulation of do-interventions #####
> doX=c(0.2, NA, NA)
> dx0.2 = dgp(10000, doX=doX, seed=SEED)
> dx0.2$df_orig$numpy()[1:5,]
     [,1]        [,2]      [,3]
[1,]  0.2 -0.06350804  1.245475
[2,]  0.2  0.23617326  4.848745
[3,]  0.2 -0.31859812 -3.545032
[4,]  0.2 -0.24405755 -6.081490
[5,]  0.2  0.53235835  7.512290
> 
> 
> doX=c(0.7, NA, NA)
> dx7 = dgp(10000, doX=doX, seed=SEED)
> #hist(dx0.2$df_orig$numpy()[,2], freq=FALSE,100)
> mean(dx7$df_orig$numpy()[,2]) - mean(dx0.2$df_orig$numpy()[,2])  
[1] -0.2000643
> mean(dx7$df_orig$numpy()[,3]) - mean(dx0.2$df_orig$numpy()[,3])  
[1] -1.170931
> 
> ########### Do(x1) seems to work#####
> 
> #### Check intervention distribution after do(X1=0.2)
> df = data.frame(train$df_orig$numpy())
> fit.x2 = Colr(X2~X1,df)
> x2_dense = predict(fit.x2, newdata = data.frame(X1 = 0.2), type = 'density')
> x2s = as.numeric(rownames(x2_dense))
> 
> ## samples from x2 under do(x1=0.2) via simulate
> ddd = as.numeric(unlist(simulate(fit.x2, newdata = data.frame(X1 = 0.2), nsim = 1000)))
> s2_colr = rep(NA, length(ddd))
> for (i in 1:length(ddd)){
+   s2_colr[i] = as.numeric(ddd[[i]]) #<--TODO somethimes 
+ }
> 
> if(sum(is.na(s2_colr)) > 0){
+   stop("Pechgehabt mit Colr, viel Glück und nochmals!")
+ }
> 
> hist(dx0.2$df_orig$numpy()[,2], freq=FALSE, 100, main='Do(X1=0.2) X2',  
+      sub='Histogram from DGP with do. Blue: Colr', xlab='samples')
> lines(x2s, x2_dense, type = 'l', col='blue', lw=2)
> 
> # fit.x3 = Colr(X3 ~ X1 + X2,df)
> # newdata = data.frame(
> #     X1 = rep(0.2, length(s2_colr)), 
> #     X2 = s2_colr)
> # 
> # s3_colr = rep(NA, nrow(newdata))
> # for (i in 1:nrow(newdata)){
> #   # i = 2
> #   s3_colr[i] = simulate(fit.x3, newdata = newdata[i,], nsim = 1)
> # }
> 
> s_dag = do_dag_struct(param_model, train$A, doX=c(0.2, NA, NA))
sample_from_target Fraction of extrapolated samples > 1 : %f 
0.217210486531258sample_from_target Fraction of extrapolated samples > 1 : %f 
0.0662188082933426> hist(dx0.2$df_orig$numpy()[,2], freq=FALSE, 50, main='X2 | Do(X1=0.2)', xlab='samples', 
+      sub='Histogram from DGP with do. red:TRAM_DAG')
> sample_dag_0.2 = s_dag[,2]$numpy()
> lines(density(sample_dag_0.2), col='red', lw=2)
> m_x2_do_x10.2 = median(sample_dag_0.2)
> 
> 
> s_dag = do_dag_struct(param_model, train$A, doX=c(0.2, NA, NA))
sample_from_target Fraction of extrapolated samples > 1 : %f 
0.214331418275833sample_from_target Fraction of extrapolated samples > 1 : %f 
0.0626999363303185> hist(dx0.2$df_orig$numpy()[,3], freq=FALSE, 50, main='X3 | Do(X1=0.2)', xlab='samples', 
+      sub='Histogram from DGP with do. red:TRAM_DAG')
> sample_dag_0.2 = s_dag[,3]$numpy()
> lines(density(sample_dag_0.2), col='red', lw=2)
> 
> 
> ###### Comparison of estimated f(x2) vs TRUE f(x2) #######
> shift_12 = shift_23 = shift1 = cs_23 = xs = seq(-1,1,length.out=41)
> idx0 = which(xs == 0) #Index of 0 xs needs to be odd
> for (i in 1:length(xs)){
+   #i = 1
+   x = xs[i]
+   # Varying x1
+   X = tf$constant(c(x, 0.5, 3), shape=c(1L,3L)) 
+   shift1[i] =   param_model(X)[1,3,2]$numpy() #2=LS Term X1->X3
+   shift_12[i] = param_model(X)[1,2,2]$numpy() #2=LS Term X1->X2
+   
+   #Varying x2
+   X = tf$constant(c(0.5, x, 3), shape=c(1L,3L)) 
+   cs_23[i] = param_model(X)[1,3,1]$numpy() #1=CS Term
+   shift_23[i] = param_model(X)[1,3,2]$numpy() #2-LS Term X2-->X3 (Beate Notation)
+ }
> 
> if (FALSE){
+   if (MA[2,3] == 'cs' && F32 == 1){
+     # Assuming xs, cs_23, and idx0 are predefined vectors
+     # Create a data frame for the ggplot
+     df <- data.frame(x2 = xs, cs_23 = cs_23)
+     
+     # Create the ggplot
+     p <- ggplot(df, aes(x = x2, y = cs_23)) +
+       geom_line(aes(color = "Complex Shift Estimate"), size = 1) +  
+       geom_point(aes(color = "Complex Shift Estimate"), size = 1) + 
+       geom_abline(aes(color = "f"), intercept = cs_23[idx0], slope = 0.3, size = 1) +  # Black solid line for 'DGP'
+       scale_color_manual(
+         values = c("Complex Shift Estimate" = "blue", "f" = "black"),  # Set colors
+         labels = c("Complex Shift Estimate", "f(x)")  # Custom legend labels with expression for f(X_2)
+       ) +
+       labs(
+         x = expression(x[2]),  # Subscript for x_2
+         y = "~f(x)",  # Optionally leave y-axis label blank
+         color = NULL  # Removes the color legend title
+       ) +
+       theme_minimal() +
+       theme(legend.position = "none")  # Correct way to remove the legend
+     
+     # Display the plot
+     p
+   } else if (MA[2,3] == 'cs' && F32 != 1){
+     # Assuming xs, shift_23, and idx0 are predefined vectors
+     # Create a data frame for the ggplot
+     df <- data.frame(x2 = xs, 
+                      shift_23 = cs_23 + ( -cs_23[idx0] - f(0)),
+                      f = -f(xs)
+                      )
+     # Create the ggplot
+     p <- ggplot(df, aes(x = x2, y = shift_23)) +
+       #geom_line(aes(color = "Shift Estimate"), size = 1) +  # Blue line for 'Shift Estimate'
+       geom_point(aes(color = "Shift Estimate"), size = 1) +  # Blue points for 'Shift Estimate'
+       geom_line(aes(color = "f", y = f), ) +  # Black solid line for 'DGP'
+       scale_color_manual(
+         values = c("Shift Estimate" = "blue", "f" = "black"),  # Set colors
+         labels = c("Shift Estimate", "f(x)")  # Custom legend labels with expression for f(X_2)
+       ) +
+       labs(
+         x = expression(x[2]),  # Subscript for x_2
+         y = "~f(x)",  # Optionally leave y-axis label blank
+         color = NULL  # Removes the color legend title
+       ) +
+       theme_minimal() +
+       theme(legend.position = "none")  # Correct way to remove the legend
+     
+     # Display the plot
+     p
+   } else{
+     print(paste0("Unknown Model ", MA[2,3]))
+   }
+  
+   
+   file_name <- paste0(fn, "_f23_est.pdf")
+   # Save the plot
+   ggsave(file_name, plot = p, width = 8, height = 8)
+   file_path <- file.path("/Users/oli/Library/CloudStorage/Dropbox/Apps/Overleaf/tramdag/figures", basename(file_name))
+   ggsave(file_path, plot = p, width = 8/3, height = 8/3)
+ }
>     
> 
> par(mfrow=c(2,2))
> plot(xs, shift_12, main='LS-Term (black DGP, red Ours)', 
+      sub = 'Effect of x1 on x2',
+      xlab='x1', col='red')
> abline(0, 2)
> 
> delta_0 = shift1[idx0] - 0
> plot(xs, shift1 - delta_0, main='LS-Term (black DGP, red Ours)', 
+      sub = paste0('Effect of x1 on x3, delta_0 ', round(delta_0,2)),
+      xlab='x1', col='red')
> abline(0, -.2)
> 
> 
> if (F32 == 1){ #Linear DGP
+   if (MA[2,3] == 'ls'){
+     delta_0 = shift_23[idx0] - f(0)
+     plot(xs, shift_23 - delta_0, main='LS-Term (black DGP, red Ours)', 
+          sub = paste0('Effect of x2 on x3, delta_0 ', round(delta_0,2)),
+          xlab='x2', col='red')
+     #abline(shift_23[length(shift_23)/2], -0.3)
+     abline(0, 0.3)
+   } 
+   if (MA[2,3] == 'cs'){
+     plot(xs, cs_23, main='CS-Term (black DGP, red Ours)', xlab='x2',  
+          sub = 'Effect of x2 on x3',col='red')
+     
+     abline(cs_23[idx0], 0.3)  
+   }
+ } else{ #Non-Linear DGP
+   if (MA[2,3] == 'ls'){
+     delta_0 = shift_23[idx0] + f(0)
+     plot(xs, shift_23 - delta_0, main='LS-Term (black DGP, red Ours)', 
+          sub = paste0('Effect of x2 on x3, delta_0 ', round(delta_0,2)),
+          xlab='x2', col='red')
+     lines(xs, -f(xs))
+   } else if (MA[2,3] == 'cs'){
+     plot(xs, cs_23 + ( -cs_23[idx0] - f(0) ),
+          ylab='CS',
+          main='CS-Term (black DGP f2(x), red Ours)', xlab='x2',  
+          sub = 'Effect of x2 on x3',col='red')
+     lines(xs, -f(xs))
+   } else{
+     print(paste0("Unknown Model ", MA[2,3]))
+   }
+ }
> #plot(xs,f(xs), xlab='x2', main='DGP')
> par(mfrow=c(1,1))
> 
> 
> if (TRUE){
+ ####### Compplete transformation Function #######
+ ### Copied from structured DAG Loss
+ t_i = train$df_orig
+ k_min <- k_constant(global_min)
+ k_max <- k_constant(global_max)
+ 
+ # from the last dimension of h_params the first entriy is h_cs1
+ # the second to |X|+1 are the LS
+ # the 2+|X|+1 to the end is H_I
+ h_cs <- h_params[,,1, drop = FALSE]
+ h_ls <- h_params[,,2, drop = FALSE]
+ #LS
+ h_LS = tf$squeeze(h_ls, axis=-1L)#tf$einsum('bx,bxx->bx', t_i, beta)
+ #CS
+ h_CS = tf$squeeze(h_cs, axis=-1L)
+ 
+ theta_tilde <- h_params[,,3:dim(h_params)[3], drop = FALSE]
+ theta = to_theta3(theta_tilde)
+ cont_dims = which(data_type == 'c') #1 2
+ cont_ord = which(data_type == 'o') #3
+ 
+ ### Continiuous dimensions
+ #### At least one continuous dimension exits
+ h_I = h_dag_extra(t_i[,cont_dims, drop=FALSE], theta[,cont_dims,1:len_theta,drop=FALSE], k_min[cont_dims], k_max[cont_dims]) 
+ 
+ h = h_I + h_LS[,cont_dims, drop=FALSE] + h_CS[,cont_dims, drop=FALSE]
+ 
+ ####### DGP Transformations #######
+ X_1 = t_i[,1]$numpy()
+ X_2 = t_i[,2]$numpy()
+ X_3 = t_i[,3]$numpy()
+ 
+ #h2 = x_2_dash = 5 * x_2 + 2 * X_1
+ h2_DGP = 5 *X_2 + 2 * X_1
+ h2_DGP_LS = 2 * X_1
+ h2_DGP_CS = rep(0, length(X_2))
+ h2_DGP_I = 5 * X_2
+ 
+ #h(x3|x1,x2) = 0.63*x3 - 0.2*x1 - f(x2)
+ h3_DGP = 0.63*X_3 - 0.2*X_1 - f(X_2)
+ h3_DGP_LS = -0.2*X_1
+ h3_DGP_CS = -f(X_2)
+ h3_DGP_I = 0.63*X_3
+ 
+ 
+ par(mfrow=c(2,2))
+ plot(h2_DGP, h[,2]$numpy(), main='h2')
+ abline(0,1,col='red')
+ confint(lm(h[,2]$numpy() ~ h2_DGP))
+ 
+ #Same for Intercept
+ plot(h2_DGP_I, h_I[,2]$numpy(), main='h2_I')
+ abline(0,1,col='red')
+ confint(lm(h_I[,2]$numpy() ~ h2_DGP_I))
+ 
+ plot(h2_DGP_LS, h_LS[,2]$numpy(), main='h2_LS')
+ abline(0,1,col='red')
+ confint(lm(h_LS[,2]$numpy() ~ h2_DGP_LS))
+ 
+ #Same for CS
+ plot(h2_DGP_CS, h_CS[,2]$numpy(), main='h2_CS')
+ abline(0,1,col='red')
+ confint(lm(h_CS[,2]$numpy() ~ h2_DGP_CS))
+ 
+ par(mfrow=c(1,1))
+ 
+ 
+ par(mfrow=c(2,2))
+ 
+ plot(h3_DGP, h[,3]$numpy(), main='h3')
+ abline(0,1,col='red')
+ confint(lm(h[,3]$numpy() ~ h3_DGP))
+ 
+ plot(h3_DGP_I, h_I[,3]$numpy(), main='h3_I')
+ abline(0,1,col='red')
+ confint(lm(h_I[,3]$numpy() ~ h3_DGP_I))
+ 
+ #same for ls  
+ plot(h3_DGP_LS, h_LS[,3]$numpy(), main='h3_LS')
+ abline(0,1,col='red')
+ confint(lm(h_LS[,3]$numpy() ~ h3_DGP_LS))
+ 
+ #same for CS
+ plot(h3_DGP_CS, h_CS[,3]$numpy(), main='h3_CS')
+ abline(0,1,col='red')
+ confint(lm(h_CS[,3]$numpy() ~ h3_DGP_CS))
+ 
+ par(mfrow=c(1,1))
+ 
+ }
> 
> 
> 
> 
> 
> 
> 
> proc.time()
  usuário   sistema decorrido 
   80.313     9.852    82.647