CODEWITHSUNDEEP
pythonarraysscipystatistics
Mike5298
Mike5298

Reputation: 385

Why does Scipy rv_discrete always return 1

I'm trying to form the distribution of a random variable. Where the probability that a reaction fires is:

Prob(Rj fires) = aj/a0

aj = a propensity function, that describes the probability an individual reaction fires in a given time period.

a0 = the sum of all the propensity functions for all the reactions in the system.

At the moment I have a function that returns an array of propensity functions for each reaction in the systems, which outputs: 

[32.  0.  0.]

I have then used a0 = sum(propensity) to calculate a0. which is followed by a second function: 

def prob_rxn_fires(propensity, a0):
    prob = propensity/a0   
    return prob

To calculate aj/a0, the output of this function is assigned to the variable: rxn_probability. There are three reactions in my system and I'm now trying to use scipy.stats.rv_discrete to sample from the distribution.

I have the number of reactions in the system stored in the array: 

num_rxn = np.arange(1, rxn_probability.size + 1).reshape(rxn_probability.shape)

which is shaped and sized so it always matches with the rxn_probability array.

My code at the moment is: 

j = stats.rv_discrete(name="Reaction index", values=(num_rxn, rxn_probability)).rvs()

but the output is always 1 and I'm not sure if thats right and if not why its not right ?

Cheers

Upvotes: 0

Views: 84

Answers (1)

FBruzzesi
FBruzzesi

Reputation: 6485

IIUC you have the following arrays:

import numpy as np

rxn_probability = np.array([1, 0, 0])

num_rxn = np.arange(1, rxn_probability.size + 1).reshape(rxn_probability.shape)
num_rxn 
array([1, 2, 3])

Therefore when using rv_discrete with values=(num_rxn, rxn_probability) you are basically saying that with probablity one it should pick the 0-th element of num_rxn.

Now if you try different values, let's say values=(num_rxn, [.8, .1, .1]) the result would be:

  • 1 with probability .8
  • 2 with probability .1
  • 3 with probability .1

If you want a uniform distribution, you can specify values=(num_rxn, [1/3] * num_rxn.size).

Upvotes: 2

Related Questions