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.