Fortran code snippets

95

John D. Cook’s math blog can often inspire a Fortran program.

Chebyshev’s inequality says that the probability of a random variable being more than k standard deviations away from its mean is less than 1/k 2. In symbols,

P(|X - E(X)| eq kigma) eq rac{1}{k^2}

This inequality is very general, but also very weak. It assumes very little about the random variable X but it also gives a loose bound. If we assume slightly more, namely that X has a unimodal distribution, then we have a tighter bound, the Vysochanskiĭ-Petunin inequality.

P(|X - E(X)| eq kigma) eq rac{4}{9k^2}

However, the Vysochanskiĭ-Petunin inequality does require that k be larger than √(8/3). In exchange for the assumption of unimodality and the restriction on k we get to reduce our upper bound by more than half.

program verify_vp
  !---------------------------------------------------------------------
  ! Verify the Vysochanskiĭ–Petunin inequality for the standard normal,
  ! and compute corresponding tail probabilities for the standardized
  ! uniform and double‐exponential (Laplace) distributions.
  !
  !   P(|Z| ≥ k) ≤ 4/(9 k^2)   for Z ~ N(0,1)
  !
  ! Tail probabilities for
  !   U ~ Uniform(–√3,√3)  (Var=1):   P(|U| ≥ k) = max(0, 1 – k/√3)
  !   L ~ Laplace(0,b) with b=1/√2 (Var=1): P(|L| ≥ k) = exp(−k/b)
  !
  !---------------------------------------------------------------------
  implicit none
  ! parameters controlling k loop
  integer, parameter :: dp = kind(1.0d0), nk = 3
  real(kind=dp), parameter :: k_min = 2.0_dp, k_inc = 1.0_dp

  ! parameters for the other distributions
  real(kind=dp), parameter :: a = sqrt(3.0_dp)         ! half‐width of Uniform(–√3,√3)
  real(kind=dp), parameter :: b = 1.0_dp/sqrt(2.0_dp)  ! scale of Laplace for Var=1
  real(kind=dp) :: k, phi_val, p_norm, bound, p_unif, p_lap
  integer :: i

  print "(A6,*(2X, A12))", "k", "Bound", " P_uniform", " P_normal", " P_laplace"

  do i=1,nk
    k = k_min + (i-1)*k_inc
    ! --- normal(0,1) tail and VP bound ---
    phi_val = 0.5_dp*(1.0_dp + erf(k/sqrt(2.0_dp)))
    p_norm  = 2.0_dp*(1.0_dp - phi_val)
    bound   = 4.0_dp/(9.0_dp*k**2)

    ! --- uniform(–√3,√3), Var=1 ---
    p_unif = merge(1.0_dp - k/a, 0.0_dp, k < a)

    ! --- Laplace(0,b), Var=1 ---
    p_lap = exp(-k/b)

    print "(F6.3,2X,*(F12.6,2X))", k, bound, p_unif, p_norm, p_lap
  end do

end program verify_vp

output:

     k         Bound     P_uniform      P_normal     P_laplace
 2.000      0.111111      0.000000      0.045500      0.059106
 3.000      0.049383      0.000000      0.002700      0.014370
 4.000      0.027778      0.000000      0.000063      0.003493

The results match Cook’s for the three distributions implemented.

2 Likes