Skip to content

Instantly share code, notes, and snippets.

@bsa7
Last active March 1, 2023 18:36
Show Gist options
  • Save bsa7/bfad22fcebce8d9310c2845c519fea7b to your computer and use it in GitHub Desktop.
Save bsa7/bfad22fcebce8d9310c2845c519fea7b to your computer and use it in GitHub Desktop.
Local and integral Laplace theorems.
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)=}')
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)}')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment