added my R scripts

1 files changed, 30 insertions, 0 deletions

diff --git a/R Scripts/likelihood.R b/R Scripts/likelihood.R
new file mode 100644
index 0000000..a4e5254
--- /dev/null
+++ b/R Scripts/likelihood.R
@@ -0,0 +1,30 @@
+# run after find-cascades
+
+load('Results/dag_dat_all.RData')
+
+edges_all = dag_dat_all
+edges_all$t2[is.na(edges_all$t2)] = max(hyp_lcc_verts$vic.day,na.rm=T)+1
+realized_all = as.numeric(rownames(dag_dat_vics)[realized])
+
+p_t = temporal(edges_all$t1, edges_all$t2, alpha)
+p_s = structural(delta, edges_all$dist)
+p = p_s * p_t
+p_tilde = 1 - p_s + p_s * temporal(edges_all$t1, edges_all$t2+1, alpha)
+
+vic.days = hyp_lcc_verts$vic.day
+vic.days[is.na(vic.days)] = max(hyp_lcc_verts$vic.day,na.rm=T)+1
+ages = vic.days-hyp_lcc_verts$spawn.date
+
+# realized edges
+prob_realized = sum(log(p[realized_all]))
+
+# non-realized edges
+prob_not_realized = sum(log(p_tilde[-realized_all]))
+
+# seeds
+prob_beta = n.seeds * log(beta)
+
+# non-seeds
+prob_beta_tilde = sum(ages) * log(1-beta)
+
+ll = prob_beta + prob_beta_tilde + prob_realized + prob_not_realized