Is there any easy way to fit probability mass function to a given dataset?
I finally figured out myself

random_array = [0,0,1,2,3]


unique, counts = np.unique(random_array, return_counts=True)


random_variable = sp.stats.rv_discrete(a = 0, b = np.inf, values = (unique, counts/np.sum(counts)))
Assume we have a list of numbers (samples)

data = [0,0,1,2,3]

I would like to fit a probability mass function for this dataset, in such a way that if I do something like

pmf.fit(data)

and by executing something like

pmf.eval(0)

I get

0.2

as return
and
by executing

pmf.eval(-1)

I get

0

as return.
Note that I am working with a discrete random variable here, so I am not fitting a pdf...