Speed up weighted average for Leela Chess Zero 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():
    nn=10000
    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.
Possible improvements:


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.
Only take variables that are 2 standard deviations from the maximum
In Monte-Carlo you spend 83% time in generating normally distributed random variables. Check lagged Fibonacci generator. Link
Monte-Carlo numpy calculations could all be combined into one loop in a lower level language (numba, C++, C).
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()


Output:

Timer unit: 3.52617e-07 s

Total time: 0.0399899 s
File: &lt;ipython-input-2-5df66cc29248&gt;
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