Local and integral Laplace theorems.

Integral.py

import math
import numpy as np

# P(k1, k2) = F(x'') - F(x')
# where x', x'' - first and second derivatives from x
# xs = (k1 - np) / (n * p * q) ** 0.5
# xss = (k2 - np) / (n * p * q) ** 0.5
# F(x) = (1 / (2 * PI) ** 0.5) * integral(0, x) e ** -(z**2 / 2) dz

def P(n, p, k1, k2):
  # Calculation of the probability that, in n trials, an event whose probability of occurrence
  # is equal to p will occur at least k1 and at most k2 times
  xs = x(k1, n, p)
  xss = x(k2, n, p)

  return F(xss) - F(xs)

def x(k, n, p):
  # Calculation x', x''
  return (k - n * p) / (n * p * (1 - p)) ** 0.5

def integral(a, b, f):
  # Calculation of the definite integral over the limits from a to b for the function f(z)
  delta = (b - a) / 100000
  result = 0
  for z in np.arange(a, b + delta, delta):
    result += f(z) * delta
    z += delta

  return result

def F(xval):
  # Finding the value of the Laplace function for arbitrary xval values
  if xval == 0:
    return 0

  def f(z):
    return math.e ** (-1 * z ** 2 / 2)

  return (1 / (2 * math.pi) ** 0.5) * integral(0, xval, f)

def FInv(t):
  # Finding the inverse Laplace function Ф(x) = t
  delta = 0.0001
  x_min = 0.0
  y_min = F(x_min)
  x_max = 10.0
  y_max = F(x_max)
  while True:
    x = (x_min + x_max) / 2
    y = F(x)

    print(f'{x_min=}, {y_min=}, {x=}, {y=}, {x_max=}, {y_max=}')
    if abs(y - t) < delta:
      break

    if y > t:
      x_max = x
    else:
      x_min = x

  return x


# Example.
# The probability of an event occurring in each of 2100 independent trials is 0.7.
# Find the probability that the event will occur: 
#   a) at least 1470 times and at most 1500 times;
#   b) at least 1470 times;
#   c) at most 1469 times

print(f'a: {P(2100, 0.7, 1470, 1500)=}')
print(f'b: {P(2100, 0.7, 1470, 2100)=}')
print(f'c: {P(2100, 0.7, 0, 1469)=}')
print(f'{FInv(0.475)=}')

Local.py

import math

# The probability that in n independent trials, in each of which
# the probability of occurrence of the event is equal to p (0 < p < 1), the event will occur
# exactly k times, no matter in what order, approximately equal to
# (the more accurate the larger n)

# Pn(k) = 1 / (n * p * q) ** 0.5 * fi(x)
# where fi(x) = 1 / (2 * PI) ** 0.5 * e ** (- 1 * x ** 2 / 2)
#     x = (k - n * p) / (n * p * q) ** 0.5

def P(n, k, p):
  # Probability calculation by local Laplace theorem
  x = (k - n * p) / (n * p * (1 - p)) ** 0.5
  return 1 / (n * p * (1 - p)) ** 0.5 * fi(x)

def fi(x):
  # Calculating phi from x
  return 1 / (2 * math.pi) ** 0.5 * math.e ** (-1 * x ** 2 / 2)

# Example
# Find the probability that event A occurs exactly 70 times in 243 trials,
# if the probability of this event occurring in each trial is 0.25

print(f'{P(243, 70, 0.25)}')