moransi.nlogo

基于netlogo仿真平台的

globals [

aaverage
moransI
Z
seed
quit
expected
serial#

]

turtles-own [

a

]

to go

  clear-all 
  ifelse fixed-seed 
    [ set seed input-seed
      random-seed seed ] 
    [ set seed new-seed 
      random-seed seed ]
  set serial# random 9999999                          
  make-turtles
  calculate-averages
  calculate-moransI
  calculate-expected
  calculate-significance
  if debug > -1 [ output ]
  set quit 1

end

to make-turtles

  ask patches [ set pcolor 98 ]
  create-turtles #agents
  ask turtles [ setxy random-pxcor random-pycor ]
  if value-distribution = "binary" [
    ask turtles [ set a random 2 ]
    ask turtles [
      if a = 0 [ set color white ]
      if a = 1 [ set color black ]
    ]
   ]
  if value-distribution = "three" [
    ask turtles [ set a random 3 ]
    ask turtles [
      if a = 0 [ set color white ]
      if a = 1 [ set color grey ]
      if a = 2 [ set color black ]
     ]
    ]    
  if value-distribution = "continuous" [
    ask turtles [ set a random-float 1 ]
    ask turtles [ set color scale-color grey a 0 1 ]
   ] 
  ask turtles [ set label-color 98 ]
  ask turtles [ set label who ]
  
end

to calculate-averages

  let atotal 0
  ask turtles [
    set atotal ( atotal + a )
  ] 
  set aaverage ( atotal / #agents ) 
  
end
  
to calculate-moransI
  
  let spatial-weight 0
  let total-spatial-weight 0
  let second-numerator-a 0
  let second-denominator-a 0

  let i 0
  let j 0

  ; adding up denominator of first term (EE wij)

  set i 0
  set j 0
  
  while [ i < #agents ] [
    set j 0                                                                   
    while [ j < #agents ] [
      
      if debug = 2 [ type "i:" type i type " j:" type j type " " type total-spatial-weight type " + " ] 
      
      set spatial-weight wij i j   
      set total-spatial-weight total-spatial-weight + spatial-weight
      
      if debug = 2 [ type spatial-weight type " = " type total-spatial-weight print ""]
      set j j + 1     
     ]
    set i i + 1
  ]
  
  if debug = 2 [ print "" ]
    
  if debug = 1 [ type "EE wij = " type total-spatial-weight print "" ]
  
  ; adding up numerator of second term ((EE wij(xi - xmean)(xj-xmean))
  
  set i 0
  set j 0 
     
  while [ i < #agents ] [
    set j 0                                                                 
    while [ j < #agents ] [
    
      ask turtle i [ set spatial-weight wij i j ]
      let i-offaverage-a ( [a] of turtle i - aaverage )
      let j-offaverage-a ( [a] of turtle j - aaverage )
      
      let ij-second-numerator-a ( spatial-weight * i-offaverage-a * j-offaverage-a )
      
      set second-numerator-a second-numerator-a + ij-second-numerator-a
      
      set j j + 1
    ]
    set i i + 1
   ]
  
  if debug = 1 [ type "EE wij(xi - xmean)(xj-xmean) = " type second-numerator-a print "" ]
  
  ; adding up denominator of 2nd term. Duplicates part of above, but simpler this way.
   
  ask turtles [                                                                   
    
    let offaverage-a ( a - aaverage )   
    set second-denominator-a second-denominator-a + ( offaverage-a * offaverage-a )
    
  ]

  if debug = 1 [ type "E(xi-xbar)^2 = " type second-denominator-a print "" print "" ]
  
  if second-denominator-a = 0 [ set quit 1 ]
  
  set moransI ( #agents / total-spatial-weight) * ( second-numerator-a / second-denominator-a )
  
end

to calculate-expected

  if formulation = "oregon"        [ set expected ( 1 / (#agents - 1) ) ]
  if formulation = "santabarbara"  [ set expected ( -1 / ( #agents - 1) ) ] 
  if formulation = "ottawa-normal" [ set expected ( -1 / ( #agents - 1) ) ] 
  if formulation = "ottawa-random" [ set expected ( -1 / ( #agents - 1) ) ] 
    
end

to calculate-significance
  
; S1,2,3,4,5 and K are all subcomponents of the rather complicated process of measuring variance
; S3 and S5 are used in oregon formulation only
; K is used in santabarbara formulation only

  let S1 0
  let S2 0
  let S3 0
  let S4 0
  let S5 0
  let K 0
    
  let i 0
  let j 0
  
  ;S1  

  let total-squared-duplicate-spatial-weight 0
  
  while [ i < #agents ] [
    set j 0                                                                 
    while [ j < #agents ] [
      
      let spatial-weight 0
          
      ask turtle i [ set spatial-weight wij i j ]
      let duplicate-spatial-weight ( spatial-weight + spatial-weight )
      let squared-duplicate-spatial-weight ( duplicate-spatial-weight * duplicate-spatial-weight )
      set total-squared-duplicate-spatial-weight total-squared-duplicate-spatial-weight + squared-duplicate-spatial-weight
      
      set j j + 1
    ]
    set i i + 1
   ]
       
  set S1 ( total-squared-duplicate-spatial-weight / 2 )

  if debug = 3 [ print "" type "S1: " type S1 print "" ]

  ;S2
  
  let total-squared-duplicate-total-spatial-weight 0
  
  set i 0
  set j 0
  
  
  while [ i < #agents ] [
    set j 0
    
    let total-spatial-weight 0
                                                                     
    while [ j < #agents ] [
        
      let spatial-weight 0  
      ask turtle i [ set spatial-weight wij i j ]
      set total-spatial-weight total-spatial-weight + spatial-weight
        
      set j j + 1
     ]
     
     let duplicate-total-spatial-weight total-spatial-weight * 2
     let squared-duplicated-total-spatial-weight duplicate-total-spatial-weight * duplicate-total-spatial-weight
     set total-squared-duplicate-total-spatial-weight total-squared-duplicate-total-spatial-weight + squared-duplicated-total-spatial-weight
          
     set i i + 1
   ]

  set S2 total-squared-duplicate-total-spatial-weight

  if debug = 3 [ type "S2: " type S2 print ""]

  ; S3  
  ; S3 is broken up into an operation on a couple of sub-summations (sigmas)
    
    let sigma1 0
    let sigma2 0
  
    ; sigma1
  
    let total-cubed-offaverage-a 0

    ask turtles [ 
      let offaverage-a a - aaverage
      let cubed-offaverage-a offaverage-a * offaverage-a * offaverage-a * offaverage-a
      set total-cubed-offaverage-a total-cubed-offaverage-a + cubed-offaverage-a
     ]
     
    set sigma1 total-cubed-offaverage-a
    if debug = 3 [ type "sigma1: " type sigma1 print ""]
    
    ; sigma2
    
    let total-squared-offaverage-a 0
    
    ask turtles [ 
      let offaverage-a a - aaverage
      let squared-offaverage-a offaverage-a * offaverage-a
      set total-squared-offaverage-a total-squared-offaverage-a + squared-offaverage-a
     ]    
    
    set sigma2 total-squared-offaverage-a
    if debug = 3 [ type "sigma2: " type sigma2 print ""]
    
    if sigma2 = 0 [ set quit 1 ]
    
  set S3 ( ( 1 / #agents ) * sigma1 ) / ( ( ( 1 / #agents ) * sigma2 ) *( ( 1 / #agents ) * sigma2 ) ) 

  if debug = 3 [ type "S3: " type S3 print "" ]

  ;S4
   
    ;sigma3
  
    let sigma3 0
    let total-spatial-weight 0  
      
    set i 0
    
    while [ i < #agents ] [
      set j 0
                                                                     
      while [ j < #agents ] [
        
        let spatial-weight 0  
        ask turtle i [ set spatial-weight wij i j ]
        set total-spatial-weight total-spatial-weight + spatial-weight
        
        if debug = 4 [ type "total-spatial-weight: " type total-spatial-weight type " + spatial-weight " type i type "->" type j type " " type spatial-weight type " = " type total-spatial-weight print "" ]
        
        set j j + 1
       ]
      
       set sigma3 total-spatial-weight 
            
       set i i + 1
     ]     
    
    if debug = 3 [ type "sigma3: " type sigma3 print ""]
  
  set S4 ( ( ( #agents * #agents ) - 3 * #agents + 3 ) * S1 - #agents * S2 + 3 * ( sigma3 * sigma3 ) )
  
  if debug = 3 [ type "S4: " type S4 print "" ]
  
  ;S5
  ; IThere is a typo in the oregon formulation, I've
  ; made a guess about how to correct it here.
  
  set S5 ( S1 - 2 * #agents * S1 + 6 * ( sigma3 * sigma3 ) )
  
  if debug = 3 [ type "S5: " type S5 print "" ]  
  
  ;K
  ;K reuses sigma1 and sigma2 from S3
    
  set K ( ( #agents * sigma1 ) / ( sigma2 * sigma2 ) )

  if debug = 3 [ type "K: " type K print "" ]  
  
  ;var
  
  let var 0
  
  if formulation = "oregon" [ 
    set var ( ( #agents * S4 - S3 * S5 ) / ( ( #agents - 1 ) * ( #agents - 2 ) * ( #agents - 3 ) * ( sigma3 * sigma3 ) ) ) ]
  
  if formulation = "santabarbara" [