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.

In [1]:
import numpy as np
from random import randint

# DONT USE CAPITALIZED NAMES FOR VARIABLES; RESERVE THESE FOR CLASSES ONLY
# ALSO, MORE DESCRIPTIVE VARIABLE NAMES WOULD HELP -- SIZE OF WHAT, FOR EXAMPLE
Size = 20
N = 75 
np.random.seed(53080)
x = Size*np.random.rand(N)
np.random.seed(23003)
y = Size*np.random.rand(N)  
#np.random.seed(10)
#Cows = np.array([randint(25,250) for p in range(N)])

# YOU ARE BETTER OFF USING NUMPY'S RANDOM NUMBER GENERATOR, SINCE YOU HAVE IMPORTED IT ALREADY
# IS THERE A REASON YOU ARE USING THE BUILTIN ONE INSTEAD?
# THEN YOU CAN SIMPLY: np.random.randint(25, 51, size=15)

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)])

# AGAIN, USING NUMPY YOU CAN CONCATENATE THESE USING np.r_ INSTEAD OF SUMMING LISTS

#np.random.seed(11)
#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)])
In [10]:
np.array([4,5,7,2,10,11]) & np.arange(6)
Out[10]:
array([0, 1, 2, 2, 0, 1])
In [7]:
#Calculates which grid square a particular location is in (turn a 2-d coordinate into a scalar)
def WhichGrid(x,y,XRange,YRange,XNum,YNum):
    #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*(XNum/XRange))*YNum+np.floor(y*(YNum/YRange))+1)
In [8]:
def Kernel(dist_squared):
    dist_squared = np.asarray(dist_squared)
    # WHY NOT JUST is_scalar = dist_squared.ndim == 0
    is_scalar = False if dist_squared.ndim > 0 else True
    # NOT CLEAR WHAT YOU ARE DOING BELOW
    dist_squared.shape = (1,)*(1-dist_squared.ndim) + dist_squared.shape
    K = 1 / (pi * (1 + dist_squared**2))
    K[(dist_squared < 0.0138)] = 0.3093
    K[(dist_squared > 60*60)] = 0
    return(K if not is_scalar else K[0])
In [9]:
from math import pi
def Iterate(Status, x, y, Suscept, Transmiss, grid, first_in_grid, last_in_grid, Num, MaxRate):
    Event = 0*Status
    INF = np.where(Status>5)[0]
    NI = INF.size # Note reported farms still infectious
    IGrids = grid[INF]-1
        
    for ii in range(NI):
        INFi = INF[ii]
        
        # THERE IS RARELY ANY NEED TO USE np.multiply. JUST USE -Transmiss[INFi]*Num
        
        trans = np.multiply(-Transmiss[INFi],Num) #transmissibility of infected farm to all other grid squares 
        maxr = MaxRate[IGrids[ii],:] #max number of animals to be infected in infected grid square
        # Elementwise multiplication
        rate = np.multiply(trans, maxr) #max number of animals to be infected in each grid square based on infected grid square
        MaxProb = 1 - np.exp(rate) #Max probability that infected farm infected noninfected farm
        rng = np.random.rand(len(MaxProb))
        m = np.where((MaxProb - rng)>0)[0]  #these grid squares need further consideration
        for n in range(len(m)):
            s = 1
            M = m[n]
            PAB = 1 - np.exp(-Transmiss[INFi]*MaxRate[IGrids[ii],M]) #Max probability that infected farm infects noninfected farms under consideration
            if (PAB == 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[INFi]*Suscept[ind1]*Kernel((x[INFi]-x[ind1])**2+(y[INFi]-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 - PAB)**(Num[M] - j)
                    if (R[j] < (PAB / 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[INFi]*Suscept[ind1]*Kernel((x[INFi]-x[ind1])**2+(y[INFi]-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 = Status + Event
    #m=np.where(Status==13); # Initiate Ring Culling Around Reported Farm
    #for i in range(len(m)):
    #    Status[m[i]]=-1;
    return {'Status':Status,'NI':NI}

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

In [10]:
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
    
    # LIBRARY IMPORT SHOULD NOT OCCUR INSIDE OF FUNCTIONS
    
    import numpy as np
    import pandas as pd
    from math import pi
    
    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
In [11]:
test[:,11]
Out[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.

In [12]:
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
In [13]:
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
Out[13]:
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

In [14]:
np.sum(long_random_everything_andind_anddist2['Status2']==1)
Out[14]:
491
In [15]:
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()
In [16]:
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
In [17]:
%matplotlib inline
import matplotlib.pyplot as plt
Matplot.summary_plot(M_pooled.intercept)
In [18]:
%matplotlib inline
import matplotlib.pyplot as plt
Matplot.summary_plot(M_pooled.first_coef)
In [19]:
%matplotlib inline
import matplotlib.pyplot as plt
Matplot.summary_plot(M_pooled.size_coef)