Problem: Given a vector of about 40 values m with normal error sd compute the weighted average of the values weighted by the chance that they are the maximum.

I have come up with 2 different approaches to solve this. The first is numerical integration, the second is monte-carlo sampling.

import numpy as np
from scipy import integrate
from scipy.stats import norm

m = np.arange(-1,1,.03)
sd = np.ones(len(m))*.12

def integ():
    n = len(m)
    dx = .05
    ts = np.arange(-2, 2, dx)
    cdfs = np.ones(len(ts))
    for i in range(n):
        cdfs *= norm.cdf(ts, m[i], sd[i])

    def func(i, x):
        return norm.pdf(x, m[i], sd[i]) / norm.cdf(x, m[i] , sd[i]) * cdfs

    ans = np.zeros(n)
    for i in range(n):
        ys = func(i, ts)
        ans[i] = integrate.simps(ys, dx=dx)
    ans *= m
    return ans

def monte():
    a = np.random.normal(size=(len(m),nn))*sd[:,None]+m[:,None]
    vals = (a==a.max(axis=0)).sum(axis=1)/nn
    vals *= m
    return vals

Both of these solutions work, but they are in the range of 1/1_000 of a second rather than the 1/1_000_000 of a second I am looking for.

This code is a prototype that I will eventually be writing in c++, but I want to make sure that this function is possible to calculate quickly enough for the change to be worth it, so I really don't care about formatting or readability, just performance.


The code below allows to calculate all 4 probabilities for top 4 values in about 110 nanoseconds in Python (using Numba).

The precision of the calculation is about +-0.02

Assumptions in the code:

  1. Standard deviation = 1 for all variables. (Different standard deviations can be easily implemented)
  2. There are always 4 values (We can assume it is top 4)

Ideas used in the code:

  1. The calculation is based on pairwise probabilities with maximum. The resulting probabilities are normalized to sum up to 1:

$$P_{i\_step1}=P(X_i>X_{max})\cdot\prod_{j\neq i,j\neq max}(1-P(X_j>X_{max}))\\P_{max\_step1}=\prod_{j\neq max}(1-P(X_j>X_{max}))\\P_i=\frac{P_{i\_step1}}{\sum_i(P_{i\_step1})}$$

  1. Normal CDF can be approximated with a sigmoid function:

$$\Phi(x) \approx \frac{1}{1+e^{-(0.07056*x^3+1.5976*x)}}$$

  1. Pairwise probabilities are calculated according to the formula:



@numba.jit(numba.float32[:,:](numba.float32[:,:]),nopython=True, nogil=True)
def prob_calc(m_ar):
    nn, n = m_ar.shape
    probs1 = np.empty(shape=(n,),dtype=np.float32)
    probs2 = np.empty(shape=(n,),dtype=np.float32)
    probs3 = np.empty(shape=(nn,n),dtype=np.float32)
    sqrt2 = numba.float32(np.sqrt(2.))
    for ii in range(nn):
        ii_maxind = np.argmax(m_ar[ii])
        ii_max = m_ar[ii,ii_maxind]
        for i in range(n):
            if i==ii_maxind:
        p_sum = numba.float32(0)
        for i in range(n):
            p = 1
            for j in range(n):
                if i == j:
                    if j!=ii_maxind:
            p_sum += p
        for i in range(n):
    return probs3

m_ar = np.random.uniform(low=-1, high=1, size=(1000,4)).astype(np.float32)

Time evaluation:

# 109 µs ± 1.48 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

This is 109 microseconds for 1000 sets of 4 values. One set is calculated in 109 nanoseconds.

I used 1000 sets to remove the influence of calling compiled code from Python. So 109 nanoseconds should all be compiled code and not the overhead time for calling a numba function from Python.

Calculated probabilities (first 5):

0.1698 0.0652 0.3825 0.3824
0.2594 0.1657 0.2844 0.2905
0.5794 0.0818 0.2350 0.1038
0.1456 0.4622 0.3049 0.0873
0.4261 0.2918 0.2210 0.0611

Real probabilities (Monte-Carlo):

0.1540 0.0411 0.4026 0.4023
0.2621 0.1481 0.2913 0.2985
0.6153 0.0633 0.2349 0.0865
0.1280 0.4929 0.3143 0.0649
0.4512 0.2965 0.2143 0.0380

Mean values:

-0.2171 -0.9985  0.4939  0.4937
 0.5553  0.1612  0.6369  0.6556
 0.8922 -0.7441  0.1040 -0.5620
-0.1951  0.8055  0.4383 -0.6080
 0.7952  0.4607  0.2183 -0.8293
  • \$\begingroup\$ Thank you so much. Do you know where most of the precision loss is? Just curious. The existence of code this fast means that this approach might work for me. \$\endgroup\$ – Oscar Smith Oct 29 '18 at 21:53
  • \$\begingroup\$ The main loss is due to probabilities not being independent. Normal cdf has a bit slower but more precise approximation via tanh. So it is not really a problem. Probabilities on the other hand are hard to predict using other means than integrating. Maybe a good ML model can do it fast, but it is a challenge. (Small neural networks with right additional features?) \$\endgroup\$ – keiv.fly Oct 29 '18 at 22:16
  • \$\begingroup\$ You can see that smaller values are in reality smaller and bigger values are bigger. So even a linear regression with a variable “difference from average” can help here. \$\endgroup\$ – keiv.fly Oct 29 '18 at 22:26

Possible improvements:

  1. Calculate for less than 40 variables. Take top 5 and calculate the probabilities for them. The Monte-Carlo shows that the probabilities go down fast. And reducing the number of variables will significantly reduce complexity.

  2. Only take variables that are 2 standard deviations from the maximum

  3. In Monte-Carlo you spend 83% time in generating normally distributed random variables. Check lagged Fibonacci generator. Link

  4. Monte-Carlo numpy calculations could all be combined into one loop in a lower level language (numba, C++, C).

  5. Use pairwise probabilities of all variables with maximum variable to approximate the probabilities.

Code for profiling Monte-Carlo:

%load_ext line_profiler
%lprun -f monte monte()


Timer unit: 3.52617e-07 s

Total time: 0.0399899 s
File: <ipython-input-2-5df66cc29248>
Function: monte at line 1

Line #      Hits         Time  Per Hit   % Time  Line Contents
     1                                           def monte():
     2         1         20.0     20.0      0.0      nn=10000
     3         1      94115.0  94115.0     83.0      a = np.random.normal(size=(len(m),nn))*sd[:,None]+m[:,None]
     4         1      19254.0  19254.0     17.0      vals = (a==a.max(axis=0)).sum(axis=1)/nn
     5         1         19.0     19.0      0.0      vals *= m
     6         1          1.0      1.0      0.0      return vals
  • \$\begingroup\$ Given how slow monte is relative to numerical integration, even if I spent 0% of my time generating normal values, this code would only be 30% faster, and I don't think I'll be able to get the time to 0 \$\endgroup\$ – Oscar Smith Oct 23 '18 at 21:29
  • \$\begingroup\$ And even if I did, that's still 500x too slow. \$\endgroup\$ – Oscar Smith Oct 23 '18 at 21:30
  • \$\begingroup\$ Won't all of these suggestions work as well with integration? \$\endgroup\$ – Oscar Smith Oct 23 '18 at 22:46
  • 1
    \$\begingroup\$ Using 1 and 2, I've got integration down to 3 seconds. I think I can do some fun math with pairwise stuff to get it down further. Thanks for your help so far. \$\endgroup\$ – Oscar Smith Oct 23 '18 at 23:28

From my performance tuning I have found that the monte-carlo solution isn't easy to make faster, but that the integration approach is. The main differences this version has it it first filters out points small enough that they will have no impact on the result. It then, computes everything except the integrals without any loops, leading to significant gains by only calling cdf and pdf once each. I could probably remove the pdf call by taking a cumulative difference of the cdf, but by now the numerical integration is taking over half the time. With these improvements, I can run this function 2500 times per second, or about 4x faster than originally. Unless anyone can speed up the integration step, I think this code is as fast as it' s going to get in python.

def integ(dx=.05):
    ind = np.argmax(m)
    m_max, sd_max = m[ind], sd[ind]
    not_small = np.where(m+4*sd > m_max)[0]
    m_big = m[not_small]
    sd_big = sd[not_small]

    n = len(m_big)
    ts = np.arange(np.min(m_big), m_max + 4*sd_max, dx)
    ts = np.transpose(np.tile(ts, (n, 1)))

    xs = (ts - m_big) / sd_big
    cdfs = norm.cdf(xs)
    prod_cdfs = np.transpose(np.tile(np.prod(cdfs, axis=1), (n, 1)))

    ys =  norm.pdf(xs) / sd_big * prod_cdfs / cdfs
    result = np.zeros(n)
    for i in range(n):
        result[i] = integrate.simps(ys[:,i], dx=dx)
    return result * m_big

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.