Journalist casualties in the Israel-Gaza war(cpj.org)

  import csv
  from collections import Counter
  from datetime import date, timedelta
  from itertools import accumulate, islice
  from matplotlib import pyplot as plt
  from scipy.stats import nchypergeom_wallenius
  from statistics import linear_regression
  
  def batched(iterable, n):
      # batched('ABCDEFG', 3) --> ABC DEF G
      if n < 1:
          raise ValueError('n must be at least one')
      it = iter(iterable)
      while batch := tuple(islice(it, n)):
          yield batch
  
  def dateab(a, b):
      # Lists dates from day a to day b.
      i = 0
      while (c := a + timedelta(days = i)) <= b:
          yield c.isoformat()
          i += 1
  
  def reader(file):
      with open(file, newline='', encoding='utf-8') as csvfile:
          for row in csv.reader(csvfile, delimiter=';', quotechar='"'):
              yield row[0:2]
  
  def sum_absdiff(iter1, iter2):
      # Compute absolute differences sum between vectors, as error measure.
      return sum(abs(x - y) for x, y in zip(iter1, iter2))
  
  days = list(dateab(date(2023, 10, 7), date(2024, 3, 1)))
  weeks = [e for e in batched(days, 7) if len(e) == 7]
  
  total_killings = [(d, int(n)) for d, n in islice(reader('gs_cummulative_total_killings.csv'), 1, None)]
  d = dict(total_killings)
  xs = list(max(d.setdefault(e, 0) for e in week) for week in weeks)
  
  journalist_killings = [t[0] for t in reader('gaza_strip_journalists_killed.csv')]
  c = Counter(journalist_killings)
  ys = list(accumulate(sum(c.setdefault(e, 0) for e in week) for week in weeks))
  
  slope, intercept = linear_regression(xs, ys, proportional=True)
  fitted = [slope * x + intercept for x in xs]
  print('Linear regression parameters: slope =', slope, ';', 'intercept =', intercept, '\n')
  
  
  M = 2375259    # 2022 population estimate (source: https://en.wikipedia.org/wiki/Gaza_Strip)
  ns = [i for i in range(1, 1200) if i % 5 == 0]
  odds = [i for i in range(1, 100) if i % 3 == 0]
  
  candidate_store = {}
  for n in ns:
      for odd in odds:
          L = []
          for x in xs:
              r = nchypergeom_wallenius(M, n, x, odd).mean()
              L.append(r)
          candidate_store[(n, odd)] = L[:]
  best_wallenius_fit = min(candidate_store, key= lambda e: sum_absdiff(candidate_store[e], ys))
  print('n, odds (best Wallenius fit):', best_wallenius_fit, '\n')
  top_20 = sorted(candidate_store, key= lambda e: sum_absdiff(candidate_store[e], ys))[:20]
  print('Top 20 (n, odds):', top_20, '\n')
  top_20_ns = [e[0] for e in top_20]
  top_20_odds = [e[1] for e in top_20]
  print('Top 20 n variation in range:', min(top_20_ns), '-', max(top_20_ns))
  print('Top 20 odds variation in range:', min(top_20_odds), '-', max(top_20_odds))
  
  fig, ax1 = plt.subplots()
  ax1.plot(xs, ys, label='Available data')
  ax1.plot(xs, fitted, label='Linear regression fit')
  ax1.plot(xs, candidate_store[best_wallenius_fit], label='Wallenius\' distribution fit')
  ax1.set_xlabel('General population (source: ochaopt.org)')
  ax1.set_ylabel('Journalists (source: wikipedia.org)')
  ax1.set_title('Cummulative killings by Israeli forces in Gaza Strip\n2023-10-07 to 2024-03-01')
  ax1.legend()
  plt.show()

  Python 3.11.8 (main, Feb  7 2024, 00:00:00) [GCC 13.2.1 20231011 (Red Hat 13.2.1-4)] on linux
  Type "help", "copyright", "credits" or "license()" for more information.
  
  ========== RESTART: XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ==========
  Linear regression parameters: slope = 0.00344597608865967 ; intercept = 0.0 
  
  n, odds (best Wallenius fit): (155, 75) 
  
  Top 20 (n, odds): [(155, 75), (160, 72), (145, 84), (165, 69), (170, 66), (150, 78), (175, 63), (140, 87), (135, 93), (150, 81), (185, 60), (140, 90), (180, 60), (190, 57), (195, 54), (200, 54), (135, 96), (180, 63), (205, 51), (185, 57)] 
  
  Top 20 n variation in range: 135 - 205
  Top 20 odds variation in range: 51 - 96