Monte Carlo Simulation of P-Value
As you mentioned, by calling np.vectorize on your monte_carlo function and applying it to your dataset b, you are essentially running a for loop over each element individually.

np.random.noncentral_chisquare takes either a float or an array of floats for its noncentrality parameter, so you could vectorize your function call by generating your entire distribution in 1 call then applying the rest of the function to the entire matrix at once.

def monte_carlo2(x, tot_sample):
    gen_dist = np.random.noncentral_chisquare(df=1, nonc=x, size=(x.shape[0],total_sample))
    compare = gen_dist &gt; x
    return np.divide(np.sum(compare, axis=1), tot_sample)

x2 = monte_carlo2(b, total_sample)


This passes your f_res data in as a vector, b, instead of 1 element at a time, generates a distribution with nonc= [b0, b1, b2, ..., bN]. This results in a gen_dist matrix of size (625527, 10) or (len(b), total_sample). b is a vertical vector, so the comparison will apply it across each row, then you just need to sum your comparison over each row (axis=1).

You could probably speed it up more by only generating one distribution per unique value of f_res, but depending on what percentage of your values are unique this may not be worth it, and I assume you want a new random sample each time.

With your original function I get a timeit of:


  15.2 s ± 125 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)


And with the vectorized version I get:


  8.84 s ± 50.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)