After comparing Keeling's Matlab code and output, and using multiple print statements throughout my code, I am still not sure why the size of the nearest infected farm is not significant. I found a couple of bugs, but they did not change the results. I also remembered that the reason I had previously decided not to use Keeling's kernel and use a Cauchy kernel was because Keeling's kernel makes the disease progress very slowly. So I have used the Cauchy kernel in the code below. The changes that I have made result in output that is consistent with Keeling's code, however I have noticed that in this new version the number of new infected is dramatically lower than it used to be.
import numpy as np
import pandas as pd
from numpy.random import randint
region_size = 20
N = 75
np.random.seed(53080)
x = region_size*np.random.rand(N)
y = region_size*np.random.rand(N)
cows = np.concatenate([randint(25, 51, size=15),
randint(51, 76, size=30),
randint(76, 95, size=20),
randint(95, 250, size=10)])
sheep = np.concatenate([randint(25, 51, size=15),
randint(51, 76, size=30),
randint(76, 95, size=20),
randint(95, 250, size=10)])
def which_grid(x, y, x_range, y_range, x_num, y_num):
'''
Calculates which grid square a particular location is in (turn a 2-d coordinate into a scalar)
Essentially: floor(Unif[0,1)griddim)griddim+floor(Unif[0,1)griddim)+1
Returns a number from 1 to griddim^2
'''
return(np.floor(x*(x_num/x_range))*y_num + np.floor(y*(y_num/y_range)) + 1)
dist_squared = 3
dist_squared = np.asarray(dist_squared)
def kernel(dist_squared):
dist_squared = np.asarray(dist_squared)
is_scalar = not dist_squared.ndim
# NOT CLEAR WHAT YOU ARE DOING BELOW
dist_squared.shape = (1,)*(1-dist_squared.ndim) + dist_squared.shape
K = 1 / (np.pi * (1 + dist_squared**2))
K[(dist_squared < 0.0138)] = 0.3093
K[(dist_squared > 60*60)] = 0
if is_scalar:
return K
return K[0]
def iterate(status, x, y, suscept, transmiss, grid, first_in_grid, last_in_grid, num, max_rate):
event = 0*status
infected = np.where(status>5)[0]
number_infected = infected.size # Note reported farms still infectious
infected_grids = grid[infected]-1
for i in range(number_infected):
infected_i = infected[i]
#transmissibility of infected farm to all other grid squares
trans = -Transmiss[infected_i] * num
#max number of animals to be infected in infected grid square
max_r = max_rate[infected_grids[i], :]
# Elementwise multiplication
# max number of animals to be infected in each grid square based on infected grid square
rate = trans * max_r
#Max probability that infected farm infected noninfected farm
max_prob = 1 - np.exp(rate)
u = np.random.rand(len(max_prob))
#these grid squares need further consideration
m = np.where((max_prob - u)>0)[0]
for n in range(len(m)):
s = 1
M = m[n]
#Max probability that infected farm infects noninfected farms under consideration
prob_infect = 1 - np.exp(-transmiss[infected_i] * max_rate[infected_grids[i], M])
if (prob_infect == 1):
# Calculate the infection probability for each farm in the susceptible grid
leng = last_in_grid[M] - first_in_grid[M] + 1
R = np.random.rand(leng)
for j in range(leng):
ind1 = first_in_grid[M]+j-1
Q = 1 - np.exp(-transmiss[infected_i]
*suscept[ind1]
*kernel((x[infected_i]-x[ind1])**2 + (y[infected_i]-y[ind1])**2))
if ((R[j] < Q) & (status[ind1] == 0)):
event[ind1] = 1
else:
R = np.random.rand(num[M])
# Loop through all susceptible farms in the grids where an infection event occurred.
for j in range(num[M]):
P = 1 - s*(1 - prob_infect)**(num[M] - j)
if (R[j] < (prob_infect / P)):
s = 0
ind1 = first_in_grid[M]+j-1
# HAVE YOU TRIED ANALYZING JUST THIS FUNCTION FOR SENSITIVITY TO
# FARM SIZE?
Q = 1 - np.exp(-transmiss[infected_i]
*suscept[ind1]
*kernel((x[infected_i]-x[ind1])**2 + (y[infected_i]-y[ind1])**2))
if ((R[j] < Q/P) & (status[ind1] == 0)):
event[ind1] = 1
# Evolve the infection process of those farms which have been exposed and already infectious ones.
status[status > 0] += 1
status += event
return {'Status':status, 'NI':number_infected}
THE FUNCTION BELOW IS CUMBERSOME. BREAK IT DOWN INTO FUNCTIONS, AND BE SURE TO WRITE TESTS TO ENSURE THAT EACH PIECE IS OPERATING AS EXPECTED
def Outbreaks(Size,N,Y0,farms,end,end2,x,y,Cows,Sheep,Maxtime=1000):
#This is an attempt of converting the Matlab Program 7.6 Code into Python
Status = np.array([0]*N) #Initial Status of each farm
init_ind = np.random.randint(0,N)
for i in range(Y0):
Status[init_ind] = 6 #one farm is initially infected
#Cows are 10.5 times more susceptible to disease than sheep
Suscept = Sheep + 10.5*Cows
Transmiss = 5.1e-7*Sheep + 7.7e-7*Cows
#Set up the grid
grid = WhichGrid(x,y,Size,Size,10.0,10.0)
tmp = sorted(grid) #Sort grid values
# DO WE REALLY NEED ALL OF THE BELOW? IF SO, BETTER TO INITIALIZE A PANDAS
# DATA FRAME THAN ALL THESE LISTS
#i = np.argsort(grid) #get indexed values after sort
i = [i[0] for i in sorted(enumerate(grid), key=lambda x:x[1])]
x = x[i]
y = y[i]
Status = Status[i]
grid = grid[i]
Transmiss = Transmiss[i]
Suscept = Suscept[i]
Cows = Cows[i]
Sheep = Sheep[i]
Xgrid = []
Ygrid = []
Num = []
first_in_grid = []
last_in_grid = []
Max_Sus_grid = []
index_inf = np.where(Status==6)[0].astype(int)
m2 = np.array(np.where(grid==1))
for i in range(1,int(max(grid))+1):
#turn the grid square number into an x-coordinate and y-coordinate (should not exceed XNum)
Xgrid.append(np.floor((i-1)/10))
Ygrid.append((i-1)%10)
m = np.array(np.where(grid==i))
Num.append(m.shape[1])
if Num[i-1] > 0:
first_in_grid.append(m.min()) #Add the "+1" here so the indicies match those in the Keeling code
last_in_grid.append(m.max())
Max_Sus_grid.append(Suscept[m].max())
else:
first_in_grid.append(0)
last_in_grid.append(-1)
Max_Sus_grid.append(0)
#Work out grid to maximum grid transmission probabilities
from numpy import ndarray
MaxRate = ndarray((max(grid),max(grid)))
#Determine maximum number of animals to be infected in each grid square
# YOUR INDENTATION IS MESSED UP HERE (only 3 spaces) WHICH SUGGESTS THIS NEVER GETS CALLED
for i in range (1,int(max(grid))):
for j in range(1,int(max(grid))):
if ((i-1)==(j-1)) | (Num[i-1]==0) | (Num[j-1] == 0):
MaxRate[i-1,j-1] = np.inf
else:
Dist2 = (Size*max([0,(abs(Xgrid[i-1]-Xgrid[j-1])-1)])/10)**2+(Size*max([0,(abs(Ygrid[i-1]-Ygrid[j-1])-1)])/10)**2
MaxRate[i-1,j-1] = Max_Sus_grid[j-1]*Kernel(Dist2)
#Susceptible, Exposed, Infectious, Reported.==> latent period is 4 days
i=1
S=len(np.where(Status==0))
E=len(np.where(np.logical_and(Status>0, Status<=5)))
I=len(np.where(np.logical_and(Status>5, Status<=9)))
R=len(np.where(Status==10))
R2=len(np.where(Status>9))
CullSheep=0
CullCattle=0
i=i+1
IterateFlag=1
# THIS IS A BAD APPROACH; GROWING LISTS IS EXTREMELY INEFFICIENT
S=[]
E=[]
I=[]
R=[]
R2=[]
CullSheep=[]
CullCattle=[]
t=[]
t.append(0)
results = np.c_[np.array([1]*N),np.arange(1,N+1),np.array([0]*N)]
# NEVER SEPARATE STATEENTS WITH SEMICOLONS; VERY HARD TO READ AND DEBUG
# JUST USE & INSTEAD OF logical_and
while(np.logical_and(t[-1]<end, IterateFlag)):
Status=Iterate(Status, x, y, Suscept, Transmiss, grid, first_in_grid, last_in_grid, Num, MaxRate)['Status']
Sus=np.where(Status==0)[0]
Exp=np.where(np.logical_and(Status>0, Status<=5))[0]
Inf=np.where(Status>5)[0]
S.append(len(Sus))
E.append(len(Exp))
I.append(len(Inf))
t.append(t[i-2]+1)
i+=1
#This is how I stop the simulation (all farms are infected)
if t[-1]>5:
if np.logical_or((E[-4]+I[-4]==0),I == N):
# JUST USE A break STATEMENT INSTEAD OF A FLAG
IterateFlag=0
# MOVE IMPORTS OUT OF FUNCTION
from scipy.stats import itemfreq
sim_num = np.array([i-1]*N)
seq = np.arange(1,N+1)
results_full = np.r_[results,np.c_[sim_num,seq,Status]]
results = results_full
#Return information regarding only farm of interest
this = results_full[np.logical_or.reduce([results_full[:,1] == x for x in farms])]
#Extract rows relating to timepoint of interest
no_this = this[this[:,0]==end]
#turn status to an indicator
Status_ind = (no_this[:,2]>5).astype(int)
# YOU SPEND A LOT OF ENERGY CONVERTING BACK AND FORTH BETWEEN NUMPY ARRAYS AND LISTS
# JUST USE ARRAYS
#Calculate distance to index farm - first infected is first in list of coords
coords = list(zip(x,y))
index = np.array((coords[index_inf][0],coords[index_inf][1]))
# THE BELOW WOULD BE MUCH CLEARER AS A LIST COMPREHENSION
dist = []
for j in range(0,N):
# CANT YOU JUST SUBINDEX THIS? coords[j,:2]
b = np.array((coords[j][0],coords[j][1]))
dist.append(np.linalg.norm(b-index))
to_return = np.c_[no_this[:,1],Status_ind,dist,Cows,Sheep,x,y]
from scipy import spatial
#Extract the infected farms
inf_farms = to_return[to_return[:,1]==1]
coords = list(zip(x,y))
#Create list of coordinates infected farms
inf_farm_coords = list(zip(inf_farms[:,5],inf_farms[:,6]))
# NOT CLEAR WHY YOU ARE DOING ALL THIS EXPENSIVE LIST CREATION
list_of_inf_coords = [list(elem) for elem in inf_farm_coords]
#Create list of coordinates of all farms
list_of_coords = [list(elem) for elem in coords]
#Calculate Euclidean distance from each farm to all infected farms- each row in matrix represents
#distance of one farm to each infected farm
dist_to_inf = spatial.distance_matrix(list_of_coords,list_of_inf_coords)
#Find distance to closest infected farm
# THERE IS A FUNCTION (np.minimum) THAT DOES ALL THIS FOR YOU
def minval(array):
#return(np.min(array[np.nonzero(array)]))
return(np.min(array))
closest_infected = np.apply_along_axis(minval,1,dist_to_inf)
average_infected = np.apply_along_axis(np.mean,1,dist_to_inf)
to_return = np.c_[to_return,closest_infected,average_infected]
#Create list of number of Cows and number of sheep for infected farms
inf_farm_cows = list(inf_farms[:,3])
inf_farm_sheep = list(inf_farms[:,4])
#Create a function that extracts farm size based on closest infected farm
def where_minval(array):
#return(np.argmin(array[np.nonzero(array)]))
return(np.argmin(array))
closest_infected_size_ind = np.apply_along_axis(where_minval,1,dist_to_inf)
closest_infected_cows = [inf_farm_cows[i] for i in closest_infected_size_ind]
closest_infected_sheep = [inf_farm_sheep[i] for i in closest_infected_size_ind]
#Returns array: farmID, Status_ind, dist_to_index, num_Cows,num_Sheep,x,y,,disttoclosestinf,avgdisttoinf,#cowsinclosestiffarm,#sheepinclosestinffarm
to_return = np.c_[to_return,closest_infected_cows,closest_infected_sheep]
#Now run the outbreak for the additional end - end2 steps
while(np.logical_and(t[-1]<end2, IterateFlag)):
Status=Iterate(Status, x, y, Suscept, Transmiss, grid, first_in_grid, last_in_grid, Num, MaxRate)['Status']
Sus=np.where(Status==0)[0]; Exp=np.where(np.logical_and(Status>0, Status<=5))[0]; Inf=np.where(Status>5)[0];
S.append(len(Sus)); E.append(len(Exp)); I.append(len(Inf));
t.append(t[i-2]+1);i+=1;
#Return only the data where a farm was not infected by the "end" day
not_inf = to_return[to_return[:,1]==0]
#append the status of these farms after the additional end2 - end days
newstatus = [Status[i] for i in not_inf[:,0]-1]
newstatus_ind = np.array([i>5 for i in newstatus]).astype(int)
final = np.c_[not_inf,newstatus_ind]
return(final)
test = Outbreaks(Size=Size,N=N,Y0=1,farms = np.arange(1,N+1),end=10,end2=50,x=x,y=y,Cows=Cows,Sheep=Sheep)
/Users/sandyalakkur/anaconda/envs/py3k/lib/python3.4/site-packages/IPython/kernel/__main__.py:56: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future /Users/sandyalakkur/anaconda/envs/py3k/lib/python3.4/site-packages/IPython/kernel/__main__.py:11: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future /Users/sandyalakkur/anaconda/envs/py3k/lib/python3.4/site-packages/IPython/kernel/__main__.py:20: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future /Users/sandyalakkur/anaconda/envs/py3k/lib/python3.4/site-packages/IPython/kernel/__main__.py:11: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future /Users/sandyalakkur/anaconda/envs/py3k/lib/python3.4/site-packages/IPython/kernel/__main__.py:20: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future /Users/sandyalakkur/anaconda/envs/py3k/lib/python3.4/site-packages/IPython/kernel/__main__.py:108: DeprecationWarning: converting an array with ndim > 0 to an index will result in an error in the future /Users/sandyalakkur/anaconda/envs/py3k/lib/python3.4/site-packages/IPython/kernel/__main__.py:160: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
test[:,11]
array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
HOW IS THE WORK BELOW DIFFERENT FROM WHAT YOU HAVE DONE ABOVE? YOU HAD ALREADY SIMULATED COWS AND SHEEP NEAR THE TOP.
import numpy as np
from random import randint
Outbreak = Outbreaks(Size=Size,N=N,Y0=1,farms = np.arange(1,N+1),end=10,end2=50,x=x,y=y,Cows=Cows,Sheep=Sheep)
Outbreak = np.c_[np.array([1]*Outbreak.shape[0]),Outbreak]
Num_outbreaks = 1000
for i in range(Num_outbreaks):
#Cows = np.array([randint(25,250) for p in range(N)])
Cows = np.array([randint(25,51) for p in range(15)]+[randint(51,76) for p in range(30)]+[randint(76,95) for p in range(20)]+[randint(95,250) for p in range(10)])
#Sheep = np.array([randint(25,250) for p in range(N)])
Sheep = np.array([randint(25,51) for p in range(15)]+[randint(51,76) for p in range(30)]+[randint(76,95) for p in range(20)]+[randint(95,250) for p in range(10)])
# ADD IS NOT A GOOD VARIABLE NAME
add = Outbreaks(Size=Size,N=N,Y0=1,farms = np.arange(1,N+1),end=10,end2=50,x=x,y=y,Cows=Cows,Sheep=Sheep)
add = np.c_[np.array([i+2]*add.shape[0]),add]
new_Outbreak = np.r_[Outbreak,add]
Outbreak = new_Outbreak
print(i,np.sum(add[:,12]==1))
/Users/sandyalakkur/anaconda/envs/py3k/lib/python3.4/site-packages/IPython/kernel/__main__.py:56: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future /Users/sandyalakkur/anaconda/envs/py3k/lib/python3.4/site-packages/IPython/kernel/__main__.py:11: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
0 3 1 0 2 0 3 0 4 0 5 0 6 0 7 1 8 0 9 0 10 1 11 1 12 0 13 1 14 1 15 0 16 1 17 1 18 0 19 0 20 1 21 0 22 2 23 0 24 0 25 1 26 0 27 1 28 1 29 0 30 0 31 0 32 0 33 3 34 3 35 0 36 0 37 0 38 0 39 0 40 2 41 2 42 0 43 0 44 0 45 0 46 0 47 0 48 0 49 0 50 0 51 1 52 0 53 0 54 1 55 2 56 3 57 1 58 2 59 0 60 0 61 0 62 2 63 0 64 2 65 0 66 1 67 0 68 1 69 2 70 0 71 0 72 0 73 0 74 1 75 0 76 1 77 1 78 0 79 0 80 0 81 0 82 0 83 0 84 0 85 2 86 0 87 0 88 0 89 0 90 0 91 0 92 0 93 3 94 0 95 0 96 0 97 1 98 1 99 1 100 0 101 0 102 0 103 0 104 1 105 0 106 0 107 0 108 0 109 1 110 2 111 0 112 1 113 1 114 0 115 2 116 0 117 0 118 1 119 0 120 0 121 0 122 0 123 0 124 1 125 1 126 2 127 1 128 0 129 1 130 0 131 0 132 0 133 0 134 0 135 0 136 0 137 0 138 2 139 1 140 2 141 1 142 1 143 0 144 0 145 1 146 0 147 0 148 2 149 0 150 1 151 0 152 0 153 0 154 1 155 0 156 2 157 1 158 2 159 1 160 0 161 0 162 0 163 0 164 0 165 0 166 0 167 0 168 0 169 0 170 3 171 3 172 0 173 0 174 0 175 0 176 0 177 0 178 1 179 0 180 1 181 1 182 1 183 0 184 0 185 0 186 0 187 0 188 0 189 0 190 0 191 0 192 1 193 1 194 0 195 1 196 2 197 0 198 0 199 0 200 1 201 0 202 0 203 0 204 4 205 0 206 0 207 1 208 1 209 0 210 0 211 0 212 1 213 0 214 1 215 0 216 0 217 0 218 0 219 0 220 0 221 0 222 0 223 0 224 0 225 0 226 0 227 1 228 4 229 0 230 0 231 0 232 1 233 0 234 3 235 0 236 0 237 1 238 0 239 0 240 2 241 0 242 0 243 0 244 1 245 0 246 0 247 0 248 1 249 0 250 1 251 0 252 0 253 0 254 4 255 0 256 1 257 2 258 0 259 0 260 0 261 2 262 0 263 0 264 0 265 1 266 0 267 3 268 0 269 0 270 1 271 0 272 0 273 1 274 0 275 4 276 0 277 1 278 0 279 0 280 0 281 0 282 0 283 0 284 0 285 0 286 0 287 0 288 0 289 2 290 0 291 1 292 0 293 1 294 1 295 0 296 1 297 1 298 0 299 0 300 0 301 0 302 0 303 0 304 2 305 0 306 0 307 1 308 1 309 0 310 2 311 0 312 1 313 1 314 0 315 1 316 0 317 2 318 1 319 2 320 2 321 1 322 0 323 0 324 2 325 2 326 0 327 1 328 0 329 0 330 2 331 0 332 1 333 1 334 0 335 1 336 0 337 1 338 2 339 0 340 1 341 0 342 0 343 0 344 0 345 0 346 2 347 0 348 0 349 0 350 1 351 0 352 1 353 0 354 1 355 0 356 0 357 3 358 0 359 0 360 1 361 0 362 0 363 0 364 0 365 0 366 0 367 0 368 1 369 0 370 0 371 0 372 0 373 2 374 0 375 2 376 2 377 0 378 0 379 0 380 1 381 0 382 0 383 1 384 0 385 0 386 0 387 4 388 1 389 0 390 0 391 2 392 0 393 1 394 0 395 2 396 0 397 0 398 1 399 0 400 1 401 0 402 0 403 0 404 1 405 0 406 0 407 0 408 0 409 0 410 0 411 2 412 0 413 0 414 2 415 1 416 0 417 0 418 0 419 2 420 0 421 0 422 0 423 1 424 0 425 1 426 1 427 0 428 2 429 0 430 1 431 1 432 2 433 0 434 2 435 0 436 1 437 0 438 0 439 1 440 0 441 0 442 0 443 0 444 1 445 0 446 0 447 0 448 0 449 0 450 0 451 0 452 1 453 0 454 0 455 0 456 0 457 0 458 0 459 3 460 0 461 1 462 1 463 0 464 0 465 0 466 0 467 0 468 1 469 0 470 1 471 0 472 0 473 0 474 0 475 0 476 1 477 0 478 0 479 0 480 0 481 2 482 0 483 2 484 0 485 0 486 1 487 0 488 0 489 0 490 0 491 1 492 0 493 0 494 0 495 0 496 0 497 0 498 1 499 1 500 0 501 2 502 0 503 4 504 1 505 1 506 0 507 0 508 0 509 1 510 1 511 0 512 2 513 2 514 0 515 1 516 1 517 0 518 0 519 0 520 1 521 0 522 1 523 0 524 0 525 1 526 2 527 0 528 0 529 1 530 0 531 0 532 0 533 0 534 0 535 0 536 1 537 0 538 0 539 0 540 0 541 0 542 1 543 0 544 1 545 0 546 0 547 0 548 0 549 2 550 4 551 1 552 1 553 0 554 1 555 1 556 2 557 2 558 1 559 0 560 0 561 0 562 0 563 0 564 1 565 0 566 0 567 0 568 0 569 1 570 0 571 1 572 1 573 0 574 0 575 0 576 0 577 0 578 0 579 0 580 0 581 2 582 0 583 0 584 0 585 0 586 0 587 0 588 1 589 0 590 0 591 0 592 0 593 1 594 0 595 0 596 0 597 0 598 0 599 3 600 1 601 2 602 1 603 0 604 0 605 1 606 0 607 0 608 0 609 0 610 0 611 0 612 0 613 0 614 0 615 0 616 0 617 0 618 1 619 1 620 0 621 0 622 1 623 0 624 0 625 0 626 2 627 1 628 0 629 0 630 1 631 0 632 0 633 0 634 0 635 1 636 0 637 0 638 0 639 1 640 0 641 0 642 1 643 0 644 0 645 1 646 1 647 1 648 0 649 3 650 1 651 0 652 0 653 0 654 0 655 1 656 1 657 0 658 0 659 0 660 1 661 1 662 0 663 0 664 0 665 0 666 1 667 0 668 2 669 0 670 2 671 0 672 0 673 1 674 1 675 1 676 0 677 1 678 2 679 1 680 0 681 2 682 0 683 1 684 0 685 2 686 0 687 0 688 0 689 0 690 1 691 0 692 0 693 2 694 0 695 0 696 1 697 0 698 0 699 2 700 0 701 0 702 0 703 0 704 0 705 2 706 2 707 0 708 0 709 0 710 0 711 0 712 0 713 0 714 0 715 0 716 0 717 0 718 0 719 0 720 1 721 0 722 1 723 0 724 3 725 0 726 0 727 0 728 0 729 0 730 2 731 1 732 0 733 0 734 1 735 0 736 0 737 1 738 1 739 0 740 0 741 1 742 1 743 0 744 0 745 0 746 0 747 0 748 2 749 0 750 2 751 1 752 2 753 2 754 0 755 0 756 0 757 0 758 0 759 0 760 1 761 0 762 1 763 0 764 0 765 1 766 0 767 0 768 0 769 1 770 0 771 0 772 0 773 0 774 0 775 1 776 2 777 0 778 0 779 2 780 0 781 0 782 0 783 1 784 0 785 0 786 0 787 0 788 0 789 2 790 0 791 0 792 0 793 1 794 1 795 0 796 0 797 0 798 0 799 0 800 0 801 1 802 2 803 0 804 0 805 0 806 0 807 0 808 1 809 1 810 1 811 1 812 1 813 1 814 1 815 0 816 2 817 0 818 1 819 0 820 1 821 0 822 0 823 1 824 0 825 1 826 3 827 1 828 0 829 0 830 0 831 1 832 0 833 0 834 0 835 0 836 0 837 0 838 0 839 1 840 0 841 0 842 0 843 0 844 0 845 0 846 0 847 0 848 0 849 0 850 1 851 0 852 0 853 0 854 1 855 0 856 0 857 1 858 0 859 0 860 2 861 0 862 2 863 0 864 0 865 0 866 0 867 0 868 0 869 1 870 0 871 1 872 3 873 0 874 0 875 0 876 0 877 0 878 0 879 0 880 1 881 0 882 0 883 2 884 0 885 1 886 1 887 1 888 0 889 0 890 1 891 1 892 0 893 0 894 2 895 1 896 0 897 0 898 3 899 1 900 1 901 0 902 0 903 1 904 1 905 0 906 2 907 0 908 1 909 0 910 3 911 1 912 0 913 0 914 1 915 1 916 0 917 0 918 2 919 1 920 0 921 1 922 0 923 0 924 1 925 2 926 0 927 0 928 0 929 0 930 0 931 0 932 1 933 0 934 1 935 0 936 0 937 0 938 0 939 0 940 0 941 2 942 0 943 0 944 0 945 0 946 0 947 0 948 0 949 0 950 0 951 0 952 1 953 0 954 1 955 0 956 0 957 1 958 0 959 1 960 0 961 0 962 2 963 0 964 0 965 0 966 1 967 1 968 1 969 2 970 0 971 1 972 0 973 0 974 0 975 2 976 0 977 0 978 0 979 1 980 0 981 1 982 0 983 2 984 0 985 0 986 0 987 0 988 0 989 0 990 1 991 0 992 0 993 1 994 0 995 1 996 0 997 0 998 2 999 0
/Users/sandyalakkur/anaconda/envs/py3k/lib/python3.4/site-packages/IPython/kernel/__main__.py:20: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future /Users/sandyalakkur/anaconda/envs/py3k/lib/python3.4/site-packages/IPython/kernel/__main__.py:108: DeprecationWarning: converting an array with ndim > 0 to an index will result in an error in the future /Users/sandyalakkur/anaconda/envs/py3k/lib/python3.4/site-packages/IPython/kernel/__main__.py:160: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
import pandas as pd
df = pd.DataFrame(Outbreak)
df.columns = ['run','farmID','Status','DistToIndex','NumCows','NumSheep','x-coord','y-coord','DistToNearestInfected','AvgDistToInfected','CowsNearestInfected','SheepNearestInfected','Status2']
long_random_everything_andind_anddist2 = df
long_random_everything_andind_anddist2
run | farmID | Status | DistToIndex | NumCows | NumSheep | x-coord | y-coord | DistToNearestInfected | AvgDistToInfected | CowsNearestInfected | SheepNearestInfected | Status2 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 1 | 1 | 0 | 20.757396 | 190 | 206 | 0.428463 | 1.864803 | 20.757396 | 20.757396 | 222 | 200 | 0 |
1 | 1 | 2 | 0 | 18.509917 | 91 | 87 | 0.391203 | 6.562081 | 18.509917 | 18.509917 | 222 | 200 | 0 |
2 | 1 | 3 | 0 | 17.417512 | 61 | 68 | 0.414427 | 10.122575 | 17.417512 | 17.417512 | 222 | 200 | 0 |
3 | 1 | 4 | 0 | 16.556293 | 43 | 31 | 0.996829 | 15.573354 | 16.556293 | 16.556293 | 222 | 200 | 0 |
4 | 1 | 5 | 0 | 16.035718 | 66 | 66 | 1.869069 | 17.528015 | 16.035718 | 16.035718 | 222 | 200 | 0 |
5 | 1 | 6 | 0 | 19.040339 | 27 | 36 | 2.842180 | 1.531988 | 19.040339 | 19.040339 | 222 | 200 | 0 |
6 | 1 | 7 | 0 | 19.803052 | 75 | 65 | 2.762186 | 0.461724 | 19.803052 | 19.803052 | 222 | 200 | 0 |
7 | 1 | 8 | 0 | 16.900562 | 65 | 69 | 2.449791 | 5.959629 | 16.900562 | 16.900562 | 222 | 200 | 0 |
8 | 1 | 9 | 0 | 16.450377 | 30 | 44 | 2.138491 | 7.732310 | 16.450377 | 16.450377 | 222 | 200 | 0 |
9 | 1 | 10 | 0 | 14.890135 | 69 | 64 | 2.561501 | 13.778480 | 14.890135 | 14.890135 | 222 | 200 | 0 |
10 | 1 | 11 | 0 | 13.484130 | 92 | 90 | 3.969138 | 13.955119 | 13.484130 | 13.484130 | 222 | 200 | 0 |
11 | 1 | 12 | 0 | 13.840768 | 67 | 64 | 3.767688 | 15.820503 | 13.840768 | 13.840768 | 222 | 200 | 0 |
12 | 1 | 13 | 0 | 16.091276 | 61 | 71 | 4.705106 | 3.921407 | 16.091276 | 16.091276 | 222 | 200 | 0 |
13 | 1 | 14 | 0 | 17.070804 | 63 | 76 | 4.246506 | 2.924207 | 17.070804 | 17.070804 | 222 | 200 | 0 |
14 | 1 | 15 | 0 | 15.546098 | 58 | 70 | 4.321338 | 5.419305 | 15.546098 | 15.546098 | 222 | 200 | 0 |
15 | 1 | 16 | 0 | 15.024453 | 70 | 54 | 5.452718 | 4.700431 | 15.024453 | 15.024453 | 222 | 200 | 0 |
16 | 1 | 17 | 0 | 16.036512 | 78 | 83 | 4.179607 | 4.741281 | 16.036512 | 16.036512 | 222 | 200 | 0 |
17 | 1 | 18 | 0 | 12.207746 | 53 | 74 | 5.860130 | 9.912863 | 12.207746 | 12.207746 | 222 | 200 | 0 |
18 | 1 | 19 | 0 | 13.352243 | 73 | 71 | 4.278867 | 11.560496 | 13.352243 | 13.352243 | 222 | 200 | 0 |
19 | 1 | 20 | 0 | 12.058908 | 34 | 28 | 5.472275 | 12.359426 | 12.058908 | 12.058908 | 222 | 200 | 0 |
20 | 1 | 21 | 0 | 11.943705 | 77 | 83 | 5.701048 | 15.881959 | 11.943705 | 11.943705 | 222 | 200 | 0 |
21 | 1 | 22 | 0 | 16.102311 | 61 | 56 | 7.859043 | 0.809430 | 16.102311 | 16.102311 | 222 | 200 | 0 |
22 | 1 | 23 | 0 | 14.762513 | 43 | 25 | 7.052037 | 3.264976 | 14.762513 | 14.762513 | 222 | 200 | 0 |
23 | 1 | 24 | 0 | 13.917072 | 91 | 76 | 6.062825 | 5.743816 | 13.917072 | 13.917072 | 222 | 200 | 0 |
24 | 1 | 25 | 0 | 11.940109 | 92 | 89 | 6.155299 | 9.874813 | 11.940109 | 11.940109 | 222 | 200 | 0 |
25 | 1 | 26 | 0 | 10.288324 | 179 | 210 | 7.325950 | 11.920267 | 10.288324 | 10.288324 | 222 | 200 | 0 |
26 | 1 | 27 | 0 | 14.432479 | 45 | 38 | 8.559388 | 2.374897 | 14.432479 | 14.432479 | 222 | 200 | 0 |
27 | 1 | 28 | 0 | 11.597777 | 42 | 36 | 9.859762 | 4.974922 | 11.597777 | 11.597777 | 222 | 200 | 0 |
28 | 1 | 29 | 0 | 11.520832 | 148 | 106 | 8.987554 | 5.926636 | 11.520832 | 11.520832 | 222 | 200 | 0 |
29 | 1 | 30 | 0 | 10.719772 | 83 | 92 | 8.402271 | 7.996031 | 10.719772 | 10.719772 | 222 | 200 | 0 |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
73976 | 1001 | 45 | 0 | 2.123537 | 89 | 86 | 12.698309 | 7.118964 | 2.123537 | 2.123537 | 67 | 54 | 0 |
73977 | 1001 | 46 | 0 | 3.459675 | 28 | 30 | 12.213462 | 8.499905 | 3.459675 | 3.459675 | 67 | 54 | 0 |
73978 | 1001 | 47 | 0 | 2.582251 | 131 | 209 | 13.633796 | 8.253553 | 2.582251 | 2.582251 | 67 | 54 | 0 |
73979 | 1001 | 48 | 0 | 4.631741 | 73 | 66 | 13.963750 | 10.385687 | 4.631741 | 4.631741 | 67 | 54 | 0 |
73980 | 1001 | 49 | 0 | 7.746719 | 107 | 106 | 13.311252 | 13.447489 | 7.746719 | 7.746719 | 67 | 54 | 0 |
73981 | 1001 | 50 | 0 | 2.718216 | 121 | 214 | 14.277536 | 3.051502 | 2.718216 | 2.718216 | 67 | 54 | 0 |
73982 | 1001 | 52 | 0 | 0.677720 | 83 | 85 | 14.668796 | 5.177794 | 0.677720 | 0.677720 | 67 | 54 | 0 |
73983 | 1001 | 53 | 0 | 3.235225 | 26 | 44 | 14.581574 | 8.995065 | 3.235225 | 3.235225 | 67 | 54 | 0 |
73984 | 1001 | 54 | 0 | 3.958684 | 75 | 60 | 15.396214 | 9.583559 | 3.958684 | 3.958684 | 67 | 54 | 0 |
73985 | 1001 | 55 | 0 | 4.084082 | 54 | 70 | 14.764359 | 9.830775 | 4.084082 | 4.084082 | 67 | 54 | 0 |
73986 | 1001 | 56 | 0 | 9.836703 | 64 | 54 | 15.416086 | 15.546453 | 9.836703 | 9.836703 | 67 | 54 | 0 |
73987 | 1001 | 57 | 0 | 14.066189 | 195 | 169 | 14.067645 | 19.832654 | 14.066189 | 14.066189 | 67 | 54 | 0 |
73988 | 1001 | 58 | 0 | 5.925386 | 76 | 73 | 17.187998 | 0.574325 | 5.925386 | 5.925386 | 67 | 54 | 0 |
73989 | 1001 | 59 | 0 | 3.610333 | 80 | 89 | 16.627714 | 2.978049 | 3.610333 | 3.610333 | 67 | 54 | 0 |
73990 | 1001 | 60 | 0 | 3.936865 | 83 | 84 | 16.576724 | 2.530989 | 3.936865 | 3.936865 | 67 | 54 | 0 |
73991 | 1001 | 61 | 0 | 3.837140 | 87 | 91 | 17.180531 | 8.346103 | 3.837140 | 3.837140 | 67 | 54 | 0 |
73992 | 1001 | 62 | 0 | 7.381665 | 69 | 65 | 17.041228 | 12.637765 | 7.381665 | 7.381665 | 67 | 54 | 0 |
73993 | 1001 | 63 | 0 | 8.545854 | 87 | 76 | 17.048799 | 13.873426 | 8.545854 | 8.545854 | 67 | 54 | 0 |
73994 | 1001 | 64 | 0 | 8.560077 | 245 | 108 | 17.451593 | 13.742620 | 8.560077 | 8.560077 | 67 | 54 | 0 |
73995 | 1001 | 65 | 0 | 14.054256 | 90 | 82 | 16.378905 | 19.674271 | 14.054256 | 14.054256 | 67 | 54 | 0 |
73996 | 1001 | 66 | 0 | 4.886189 | 38 | 51 | 18.330410 | 2.952671 | 4.886189 | 4.886189 | 67 | 54 | 0 |
73997 | 1001 | 67 | 0 | 5.822442 | 80 | 90 | 18.864681 | 2.107629 | 5.822442 | 5.822442 | 67 | 54 | 0 |
73998 | 1001 | 68 | 0 | 4.509184 | 63 | 54 | 18.775315 | 6.568489 | 4.509184 | 4.509184 | 67 | 54 | 0 |
73999 | 1001 | 69 | 0 | 4.066472 | 54 | 73 | 18.186413 | 7.081517 | 4.066472 | 4.066472 | 67 | 54 | 0 |
74000 | 1001 | 70 | 0 | 6.158406 | 69 | 56 | 18.492750 | 10.314425 | 6.158406 | 6.158406 | 67 | 54 | 0 |
74001 | 1001 | 71 | 0 | 6.034976 | 52 | 74 | 18.078301 | 10.504857 | 6.034976 | 6.034976 | 67 | 54 | 0 |
74002 | 1001 | 72 | 0 | 9.298684 | 72 | 69 | 19.985391 | 13.156041 | 9.298684 | 9.298684 | 67 | 54 | 0 |
74003 | 1001 | 73 | 0 | 10.615504 | 93 | 81 | 18.989877 | 15.310801 | 10.615504 | 10.615504 | 67 | 54 | 0 |
74004 | 1001 | 74 | 0 | 9.998750 | 83 | 93 | 19.700055 | 14.208172 | 9.998750 | 9.998750 | 67 | 54 | 0 |
74005 | 1001 | 75 | 0 | 12.295354 | 76 | 92 | 18.040890 | 17.493436 | 12.295354 | 12.295354 | 67 | 54 | 0 |
74006 rows × 13 columns
np.sum(long_random_everything_andind_anddist2['Status2']==1)
491
dist_diff = (long_random_everything_andind_anddist2['DistToNearestInfected'] - long_random_everything_andind_anddist2['DistToNearestInfected'].mean())/5
closest_cow_diff = (long_random_everything_andind_anddist2['CowsNearestInfected'] - long_random_everything_andind_anddist2['CowsNearestInfected'].mean())/20
closest_sheep_diff = (long_random_everything_andind_anddist2['SheepNearestInfected'] - long_random_everything_andind_anddist2['SheepNearestInfected'].mean())/20
status = long_random_everything_andind_anddist2['Status2']
from pymc import Normal, Binomial, Gamma, Lambda, invlogit, MCMC, Matplot, Bernoulli, MAP, AdaptiveMetropolis
#N = df.shape[0]
def pooled_model():
# Common slope & intercept prior
intercept = Normal('intercept', mu=0., tau=0.001, value = 0)
first_coef = Normal('first_coef', mu=0., tau=0.001, value = 0)
size_coef = Normal ('size_coef', mu=0., tau= 0.001, value=[0]*2)
#likelihood model
prob = Lambda('prob', lambda intercept=intercept,first_coef=first_coef, size_coef= size_coef:
invlogit(intercept + first_coef*dist_diff + size_coef[0]*closest_cow_diff +
size_coef[1]*closest_sheep_diff))
y = Bernoulli('y', p=prob, value=status, observed=True)
return locals()
chains = 2
iterations = 10000
burn = 4000
M_pooled = MCMC(pooled_model())
M_map = MAP(pooled_model())
M_pooled = MCMC(M_map)
#M_pooled.use_step_method(AdaptiveMetropolis, M_pooled.b)
for i in range(chains):
M_pooled.sample(iterations, burn)
[-----------------100%-----------------] 10000 of 10000 complete in 128.1 sec
%matplotlib inline
import matplotlib.pyplot as plt
Matplot.summary_plot(M_pooled.intercept)
%matplotlib inline
import matplotlib.pyplot as plt
Matplot.summary_plot(M_pooled.first_coef)
%matplotlib inline
import matplotlib.pyplot as plt
Matplot.summary_plot(M_pooled.size_coef)