sim_doctors <- function(num.doctors, num.days, initial_doctors, p) {
# Remember to set up all_doctors
all_doctors <- 1:num.doctors
# Remember to set up has_adopted as binary vector
has_adopted <- matrix(0,nrow=num.doctors,ncol=num.days)
# Set some doctors to have initially adopted
# initial_doctors are indices of doctors who are using as of day 1
has_adopted[initial_doctors,1:num.days] <- 1
for (today in 1:num.days) {
# pull two random doctors
todays_doctors <- sample(all_doctors,size=2,replace=FALSE)
# check that one has adopted and the other hasn't
if(sum(has_adopted[todays_doctors,today])==1) {
# make the non-adopter adopt with probability p
which_of_todays <- which(has_adopted[todays_doctors,today]==0)
receiver <- todays_doctors[which_of_todays]
has_adopted[receiver,today:num.days] <- rbinom(n=1,size=1,prob=p)
}
}
return(has_adopted)
}
sim1 <- sim_doctors(num.doctors=10,num.days=1000,initial_doctors=3,p=1)
plot(colSums(sim1))
sim2 <- sim_doctors(num.doctors=100,num.days=1000,initial_doctors=3,p=0.5)
plot(colSums(sim2))