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
|
library(igraph)
setwd("~/Documents/Violence Cascades/")
load('Results/hyper-lcc.RData')
# load data
hyp_lcc = upgrade_graph(hyp_lcc)
vic_ids = which(V(hyp_lcc)$vic==TRUE)
non_vic_ids = which(V(hyp_lcc)$vic==FALSE)
hist(as.numeric(hyp_lcc_verts$vic.day[vic_ids]),2000,col='lightblue',
xlab='Day of Study Period',main='Infections During the Study Period')
# get infection counts per day
days = 1:3012
t = table(factor(hyp_lcc_verts$vic.day[vic_ids],levels=days))
counts = as.vector(t)
infs = data.frame(days,counts)
# define background function
fit = function(x, mu, A, phi) {mu + A*(sin((pi/365)*x+phi))^2}
fit_form = counts ~ mu + A*(sin((pi/365)*days+phi))^2
# explore data
plot(days,counts)
curve(fit(x, mu=1, A=12, phi=2.8), add=TRUE ,lwd=2, col="steelblue")
res = nls(formula=fit_form, data=infs, start=list(mu=1, A=12, phi=2.8))
co = coef(res); co
plot(t)
plot(days,counts,pch=20,cex=0.4,xlab='Day',ylab='Number of Infections')
curve(fit(x, mu=co["mu"], A=co["A"], phi=co["phi"]), add=TRUE ,lwd=5, col="steelblue")
|
