Question: Generating bivariate zero-inflated discrete random numbers

I want to generate bivaraite zero-inflated descrite random numbers from this distribution defined below

### Generating Random Numbers from MBZIPR`  

restart:      

with(Statistics):  

randomize():  

N    := 100;

`x__01` := Vector[row](P, [1$N]);  

`x__11` := Sample(Binomial(N, 0.4), N);  

`x__22` := Sample(Normal(0, 1), N);  

`z__02` := Vector[row](P, [1$N]);  

`z__12` := Sample(Binomial(N, 0.4), N);  

`z__22` := Sample(Normal(0, 1), N);  

phi__1:= (1)/((1+exp(-(gamma__01 *`z__01` +gamma__11*`z__11` + gamma__21*`z__21`)))):      

phi__2:=(1)/((1+exp(-(gamma__02 *`z__02` +gamma__12*`z__12` + gamma__22*`z__22`)))):    

lambda__1:=exp(beta__01 *`x__01` +beta__11*`x__11` + beta__21*`x__21`)*(1+exp(gamma__01 *`z__01` +gamma__11*`z__11` + gamma__21*`z__21`)):    

lambda__2:=exp(beta__02 *`x__02` +beta__12*`x__12` + beta__22*`x__22`)*(1+exp(gamma__02 *`z__02` +gamma__12*`z__12` + gamma__22*`z__22`)):    

B:= (y[1],  y[2])->([[phi +(1-phi)*((e)^(-lambda[1]- lambda[2])*(1+ alpha*(1-(e)^(-(1-(e)^(-1))*lambda[1]))*(1-(e)^(-(1-(e)^(-1))*lambda[2])))),],[(1-phi)*((e)^(-lambda[1]- lambda[2])*((lambda[1])^(y[1]) *  (lambda[2])^(y[2]))/(y[1]!* y[2]!)*((1+ alpha)*((e)^(-y[1])-(e)^(-(1-(e)^(-1))*lambda[1]))*((e)^(-y[2])-(e)^(-(1-(e)^(-1))*lambda[2])))),]]):  

B1:= Statistics:-Distribution(  Type= discrete,        ProbabilityFunction= B,        Support= 0..infinity,        DiscreteValueMap= (n-> n)  

):  (beta__01,beta__11,beta__21,beta__02,beta__12,beta__22,gamma__01,gamma__11,gamma__21,gamma__02,gamma__12,gamma__22,phi,alpha):= (0.2,  -2,  0.25,  0.15,  -2.5,  0.2, 0.1,  2,  -2.5,  0.3,  1.3,  2.5,  0.5, -2) :  M:= Matrix((100, 25), datatype= float[8]);  S:=Statistics:-Sample(B1, M, method= [discrete, range= 0..100]);