```
import math
import random
from scipy.integrate import quad

assumed_incident_duration = 3
rate = 100

def instant_bonus(hour_start: float) -> float:
    return (8 / math.sqrt(2 * math.pi)) * math.exp((-0.5*((hour_start - 4))**2)) + 1.5

def total_bonus(time_start: float, hours_spent: float, rate: int) -> float:
    time_end = time_start + hours_spent
    bonus = rate * instant_bonus(time_start) + rate * quad(instant_bonus, time_start, time_end)[0]
    return bonus

def random_bonus(incident_duration: float) -> float:
    uniform_end = 8-incident_duration
    start_time = random.uniform(0, uniform_end)
    return total_bonus(start_time, incident_duration, rate)

sum_bonus = 0
for i in range(100000):
    sum_bonus += random_bonus(assumed_incident_duration)

print(f"average bonus = {sum_bonus / 100000}")
# average bonus = 1187.025968163705
# almost always between 1187 and 1189
```