Monte Carlo Simulation of P-Value

Question

I'm testing Python 3 code to perform a Monte Carlo simulation based on the result of an statistical test.

I currently have the result of the statistical test in a pandas dataframe, like this.

Dataframe A
+-----+-------+
| id  | f_res |
+-----+-------+
|   1 | 4.22  |
|   2 | 5.25  |
|   3 | 3.3   |
|   4 | 2.5   |
|   5 | 1.9   |
|   6 | 9.3   |
+-----+-------+

So my idea is; for each row in f_res, pass that value to a function and extract multiple values from a Noncentral chi-squared distribution, ask how many of this extracted values are greater than the original and divide this by the total of values analyzed.

I'm working with numpy to generate the array of values, and I have this working code.

import numpy as np
import pandas as pd

total_sample = 100

def monte_carlo(x, tot_sample):
    gen_dist = np.random.noncentral_chisquare(df=1, nonc=x, size=tot_sample)
    compare = gen_dist > x
    return np.divide(np.sum(compare), tot_sample)

df = pd.DataFrame(np.random.randint(0,10, 625527), columns=[['A']])
b = np.array(df.values)
g = np.vectorize(monte_carlo)

x = g(b, total_sample)

I know that the use of np.vectorize is not good for speed, because behind curtains is a for loop, but currently I don't understand another method for applying a function element wise in numpy.

Currently analyzing a dataframe of 625.000 items and with a total sample of 100, it takes approx 20 s. line_profiler shows the following.

Total time: 19.909 s File: test3.py Function: test at line 6

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     6                                           @profile
     7                                           def monte_carlo(x, tot_sample):
     8    625528     33866564     54.1     58.2      gen_dist = np.random.noncentral_chisquare(df=1, nonc=x, size=tot_sample)
     9    625528      4621263      7.4      7.9      compare = gen_dist > x
    10    625528     19704575     31.5     33.9      return np.divide(np.sum(compare), tot_sample)

So, using numpy or pandas or both, is there any way I could improve this in execution time and syntax?

The final idea is to return a numpy array with the simulated p-values to append it to the original dataframe.

Answer 1

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 > 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)