created examples folder

ucl-cssb · Nov 27, 2023 · e972169 · e972169
1 parent 4df9c3d
commit e972169
Show file tree

Hide file tree

Showing 10 changed files with 7,987 additions and 0 deletions.
diff --git a/examples/examples-Stein.ipynb b/examples/examples-Stein.ipynb
diff --git a/examples/examples-mcmc-maria.ipynb b/examples/examples-mcmc-maria.ipynb
diff --git a/examples/examples-mcmc-ode.ipynb b/examples/examples-mcmc-ode.ipynb
diff --git a/examples/examples-mcmc.ipynb b/examples/examples-mcmc.ipynb
diff --git a/examples/examples-sde-maria.ipynb b/examples/examples-sde-maria.ipynb
diff --git a/examples/examples-simulated.ipynb b/examples/examples-simulated.ipynb
diff --git a/examples/model_nopert.txt b/examples/model_nopert.txt
@@ -0,0 +1,74 @@
+functions {
+  vector lotka_volterra_N_red(real[] x, int N, vector mu, vector Md, vector M) {
+     // Model: F = dlnX/dt = mu + M x
+
+     vector[N] dydt;
+
+     int countM = 1;
+
+     for(i in 1:N){
+        dydt[i] = mu[i] - Md[i]*x[i];
+
+        for(j in 1:N){
+            if ( i != j ){
+                dydt[i] += M[countM]*x[j];
+                countM += 1; 
+                //print("loop iteration: ", i, j, countM);
+             }
+         }
+     }
+
+     return dydt;
+  }
+
+}
+
+data {
+  int<lower=1> N;
+  int<lower=1> T;
+
+  array[T,N] real y;
+  array[T,N] real x;
+
+  real tau0;
+  real sigma;
+}
+
+parameters {
+  vector<lower=0>[N]  mu;
+  vector<lower=0>[N]  Md;
+  vector[N*N - N]     M;
+
+  vector<lower=0>[N*N - N]  lambda;
+  real<lower=0>  tau;
+
+  //real<lower=0>  sigma;
+}
+
+model {
+
+  // mu
+  target += lognormal_lpdf(mu | 0.01, 0.5);
+
+  // Md 
+  target += normal_lpdf(Md | 0.1, 0.05);
+
+  // Mij: Horsehoe prior
+  target += cauchy_lpdf(tau | 0, tau0);
+
+  for(i in 1:(N*(N-1))){
+        target += normal_lpdf(M[i] | 0, lambda[i]*tau);
+        target += cauchy_lpdf(lambda[i] | 0, 1);
+  }
+  //
+
+  // sigma
+  //target += lognormal_lpdf(sigma | 0.01, 0.5);
+
+  for (t in 1:T) {
+      vector[N] y_hat = lotka_volterra_N_red(x[t,:], N, mu, Md, M);
+      for (s in 1:N){
+        target += normal_lpdf(y[t,s] | y_hat[s], sigma);
+      }
+    }
+}
diff --git a/examples/model_pert.txt b/examples/model_pert.txt
@@ -0,0 +1,85 @@
+
+functions {
+  vector lotka_volterra_N_red(real[] x, int N, vector mu, vector Md, vector M, vector E, real u) {
+     // Model: F = dlnX/dt = mu + M x + E u
+
+     vector[N] F;
+
+     int countM = 1;
+
+     for(i in 1:N){
+        F[i] = mu[i] - Md[i]*x[i];
+
+        // off diagonal interaction terms 
+        for(j in 1:N){
+            if ( i != j ){
+                F[i] += M[countM]*x[j];
+                countM += 1; 
+                //print("loop iteration: ", i, j, countM);
+             }
+         }
+
+         // epsilon terms
+         F[i] += E[i]*u;
+     }
+
+     return F;
+  }
+
+}
+
+data {
+  int<lower=1> N;
+  int<lower=1> T;
+  int<lower=0> Np;
+
+  array[T,N]  real y;
+  array[T,N]  real x;
+  array[T,Np] real u;
+
+  real sigma;
+  real tau0;
+}
+
+parameters {
+  vector<lower=0>[N]  mu;
+  vector<lower=0>[N]  Md;
+  vector[N*N - N]     M;
+  vector[N]           E;
+
+  vector<lower=0>[N*N - N]  lambda;
+  real<lower=0>  tau;
+
+  //real<lower=0>  sigma;
+}
+
+model {
+
+  // mu
+  target += lognormal_lpdf(mu | 0.01, 0.5);
+
+  // Md 
+  target += normal_lpdf(Md | 0.1, 0.05);
+
+  // Mij: Horsehoe prior
+  target += cauchy_lpdf(tau | 0, tau0);
+
+  for(i in 1:(N*(N-1))){
+        target += normal_lpdf(M[i] | 0, lambda[i]*tau);
+        target += cauchy_lpdf(lambda[i] | 0, 1);
+  }
+  //
+
+  // sigma
+  //target += lognormal_lpdf(sigma | 0.01, 0.5);
+
+  // epsilon
+  target += normal_lpdf(E | 0, 0.5);
+
+  for (t in 1:T) {
+      vector[N] y_hat = lotka_volterra_N_red(x[t,:], N, mu, Md, M, E, u[t,1] );
+      for (s in 1:N){
+        target += normal_lpdf(y[t,s] | y_hat[s], sigma);
+      }
+    }
+}
diff --git a/examples/model_sde_nopert.txt b/examples/model_sde_nopert.txt
@@ -0,0 +1,86 @@
+functions {
+  vector gLV_sde(real[] x, int N, vector mu, vector Md, vector M) {
+     // Model: f(x) = mu + M x
+
+     vector[N] f;
+
+     int countM = 1;
+
+     for(i in 1:N){
+        f[i] = mu[i] - Md[i]*x[i];
+
+        for(j in 1:N){
+            if ( i != j ){
+                f[i] += M[countM]*x[j];
+                countM += 1; 
+                //print("loop iteration: ", i, j, countM);
+             }
+         }
+     }
+
+     return f;
+  }
+
+}
+
+data {
+  int<lower=1> N;
+  int<lower=1> T;
+
+  array[T,N] real x;
+  array[T] real times;
+  array[T,N] int xmiss;
+
+  real tau0;
+  //real sigma;
+}
+
+parameters {
+  vector<lower=0>[N]  mu;
+  vector<lower=0>[N]  Md;
+  vector[N*N - N]     M;
+
+  vector<lower=0>[N*N - N]  lambda;
+  real<lower=0>  tau;
+
+  real<lower=0>  sigma;
+
+  //real d0;
+  //real d1;
+}
+
+model {
+
+  // mu
+  target += lognormal_lpdf(mu | 0.01, 0.5);
+
+  // Md 
+  target += normal_lpdf(Md | 0.1, 0.05);
+
+  // Mij: Horsehoe prior
+  target += cauchy_lpdf(tau | 0, tau0);
+
+  for(i in 1:(N*(N-1))){
+        target += normal_lpdf(M[i] | 0, lambda[i]*tau);
+        target += cauchy_lpdf(lambda[i] | 0, 1);
+  }
+
+  // Negative binomial observation error
+  //target += neg_binomial_lpmf(x[i] | x[i], d0/x[i] + d1)
+
+  // sigma
+  target += lognormal_lpdf(sigma | 0.01, 0.5);
+
+  // SDE likelihood
+  // log xi(t+1) ~ N( log xi(t) + f(x)dt,  )
+
+  for (t in 1:(T-1)) {
+      vector[N] f_x = gLV_sde(x[t,:], N, mu, Md, M);
+      real dt = times[t+1] - times[t];
+
+      for (s in 1:N){
+         target += normal_lpdf(x[t+1,s] | x[t,s] + f_x[s]*dt, sqrt(dt)*sigma);
+      }
+  }
+  // tmp
+}