From 1739e9f5706bb8a73de5dbf0b467de49ea040898 Mon Sep 17 00:00:00 2001
From: Ben Green <bgreen@g.harvard.edu>
Date: Mon, 8 Jun 2015 15:21:51 -0400
Subject: added my R scripts

---
 R Scripts/likelihood.R | 30 ++++++++++++++++++++++++++++++
 1 file changed, 30 insertions(+)
 create mode 100644 R Scripts/likelihood.R

(limited to 'R Scripts/likelihood.R')

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
-- 
cgit v1.2.3-70-g09d2