Advanced Search

Josolumoh

25 Reputation

5 Badges

9 years, 162 days
Contact Josolumoh

MaplePrimes Activity


These are replies submitted by Josolumoh

The goal is to generate bivariate random...

Josolumoh 25
November 19 2023
0 0

@mmcdara

I have effected all the suggetions. Please help with last stage.  

restart;
with(Statistics);
randomize();
N := 100;
x__01 := Vector[row](P, [1 $ N]);
x__02 := Vector[row](P, [1 $ N]);
x__11 := Sample(Binomial(N, 0.4), N);
x__12 := Sample(Binomial(N, 0.4), N);
x__21 := Sample(Normal(0, 1), N);
x__22 := Sample(Normal(0, 1), N);
z__01 := Vector[row](P, [1 $ N]);
z__02 := Vector[row](P, [1 $ N]);
z__11 := Sample(Binomial(N, 0.4), N);
z__12 := Sample(Binomial(N, 0.4), N);
z__21 := Sample(Normal(0, 1), N);
z__22 := Sample(Normal(0, 1), N);
t := 1, 2;
aux1 := 1 +~ exp~(-((`*`~(gamma__01, z__01) +~ `*`~(gamma__11, z__11)) +~ `*`~(gamma__21, z__21)));
aux2 := 1 +~ exp~(-((`*`~(gamma__02, z__02) +~ `*`~(gamma__12, z__12)) +~ `*`~(gamma__22, z__22)));
phi__1 := `/`~(1, aux1);
phi__2 := `/`~(1, aux2);
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 -> [[phi__t +~ (1 - phi__t) *~ exp~(-lambda[1] - lambda[2]) *~ (1 +~ `*`~(alpha, 1 - exp~(-`*`~(1 - exp~(-1), lambda[1]))) *~ (1 - exp~(-(1 - exp~(-1))*lambda[2])))], [(1 - phi__t) *~ exp~(-lambda[1] - lambda[2]) *~ `*`~(lambda[1]^y[1], lambda[2]^y[2])/y[1]! *~ y[2]! *~ `+`~(1, alpha) *~ (exp~(-y[1]) - exp(-`*`~(1 - exp(-1), lambda[1])))*(exp(-y[2]) - exp(-`*`~(1 - exp(-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, alpha := 0.2, -2, 0.25, 0.15, -2.5, 0.2, 0.1, 2, -2.5, 0.3, 1.3, 2.5, -2;
M := Matrix(100, 25, datatype = float[8]);
S := Statistics:-Sample(B1, M, method = [discrete, range = 0 .. 100]);

Generating Discrete Random Numbers from ...

Josolumoh 25
July 19 2023
0 0

I am aware of the robustness of the Maple in generating randoms for simulation study.

I really appreciate the example of the zero-inflated Poisson given; I think ours should follow the approach but, I don’t know I to  do it especially when it is bivariate.   

Please, I need help in generating discrete random numbers form my newly developed zero-inflated distribution defined in B:

 

#`### Generating Random Numbers from MBZIPR`  

restart:    

with(Statistics):  

randomize():  

N    := 100; 

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

x__1 := Sample(Binomial(N, 0.4), N);

x__2 := Sample(Normal(0, 1), N);  

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

z__1 := Sample(Binomial(N, 0.4), N);

z__2 := Sample(Normal(0, 1), N) ;    

t:= (1, 2);  

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__t +(1-phi__t)*(exp(-lambda[1]- lambda[2])*(1+ alpha*(1-exp(-(1-(e)^()xp(-1))*lambda[1]))*(1-exp(-(1-exp(-1))*lambda[2])))),],[(1-phi__t)*(exp(-lambda[1]- lambda[2])*((lambda[1])^(y[1]) *  (lambda[2])^(y[2]))/(y[1]!* y[2]!)*((1+ alpha)*(exp(-y[1])-exp(-(1-exp(-1))*lambda[1]))*(exp(-y[2])-exp(-(1-exp(-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]); 

Many thanks for the support.

Thanks...

Josolumoh 25
July 19 2023
0 0

@acer 

Really appreciate. 

Error, (in unapply)...

Josolumoh 25
March 17 2019
0 0


Thank you Prof,

 

I am very grateful for the supports.

 

In trying to run the other part I ran in to error 

Error, (in unapply) variables must be unique and of type name

  

 

How do I apply the b-values correctly ?

 

 

 

 

 

 

> 

restart; with(Statistics); randomize(); B := unapply*(simplify*binomial(x+r-1, x)*((1/2)*b/(d*x+1+(1/2)*b))^r*((d*x+1)/(d*x+1+(1/2)*b))^x/(d*x+1), [r, b, d]); beta := Array(0 .. 2, [.25, .6, -.2]); N := 5

155280730409912

  

unapply*(simplify*binomial(x+r-1, x)*((1/2)*b/(d*x+1+(1/2)*b))^r*((d*x+1)/(d*x+1+(1/2)*b))^x/(d*x+1), [r, b, d])

  

beta := Array(0..2, {(1) = .25, (2) = .6})

  

5

 (1)
> 

X__0 := 1; X__1 := RandomVariable(Bernoulli(1/2)); X__2 := RandomVariable(Normal(0, 1)); SX__1 := Sample(X__1, N); SX__2 := Sample(X__2, N)

`#msub(mi("X"),mi("0"))` := 1

  

`#msub(mi("X"),mi("1"))` := _R

  

`#msub(mi("X"),mi("2"))` := _R0

  

`#msub(mi("SX"),mi("1"))` := Vector[row](5, {(1) = 1.0, (2) = 1.0, (3) = .0, (4) = .0, (5) = .0}, datatype = float[8])

  

`#msub(mi("SX"),mi("2"))` := Vector[row](5, {(1) = -.6899027608083386, (2) = .4513522093757506, (3) = .4441150943533687, (4) = -.3294600013068973, (5) = -1.734357540085828}, datatype = float[8])

 (2)
> 

X__B := proc (x, y, z) options operator, arrow; exp(add(beta[k]*[x, y, z][k+1], k = 0 .. 2)) end proc;

proc (x, y, z) options operator, arrow; exp(add(beta[k]*[x, y, z][k+1], k = 0 .. 2)) end proc

 (3)
> 

Sample_b := Sample(X__B(X__0, X__1, X__2), N)

Sample_b := Vector[row](5, {(1) = 2.537110948198722, (2) = 2.2344320592161195, (3) = 2.0967402658713548, (4) = 1.7035732378545358, (5) = 2.558742697235122}, datatype = float[8])

 (4)
> 

NULL

[HFloat(2.685805148674752), HFloat(2.137698019254597), HFloat(1.1748928842842696), HFloat(1.3714821279159721), HFloat(1.8164272239994432)]

 (5)
> 

> 

NULL
> 

SUBS := proc (n) options operator, arrow; [r = 1, b = val_b[n], d = 0] end proc;

proc (n) options operator, arrow; [r = 1, b = val_b[n], d = 0] end proc

 (6)
> 

Xmax := 25:
> 

NperB := 50:
> 

U := Statistics:-RandomVariable(Uniform(0, 1)):
> 

AllSamples := NULL:
> 

for n to N do PointwiseB := [seq([x, B(`~`[rhs](SUBS(n))[])], `in`(x, [`$`(0 .. Xmax)]))]; ShouldBe1 := add(`~`[op](2, PointwiseB)); PointwiseB := [seq([x, B(`~`[rhs](SUBS(n))[])/ShouldBe1], `in`(x, [`$`(0 .. Xmax)]))]; c := CumulativeSum(`~`[op](2, PointwiseB)); ApproximateCDF := zip(proc (u, v) options operator, arrow; [u[1], v] end proc, PointwiseB, convert(c, list)); AllSamples := AllSamples, map(proc (u) options operator, arrow; ListTools:-BinaryPlace(`~`[op](2, ApproximateCDF), u) end proc, Sample(U, NperB)) end do; AllSamples := `<,>`(AllSamples)

Error, (in unapply) variables must be unique and of type name

  

 

AllSamples := Matrix(0, 1, {})

 (7)
> 

``


 

Download Second_Scenario.mw

@tomleslie 

Thank You...

Josolumoh 25
March 16 2019
0 0

@tomleslie 

Hi Prof,

 

X__B:= (x, y, z)-> exp(add(beta[k]*[x,y,z][k+1], k = 0 .. 2));

This returns the expeceted results.

I will now try and implement it in the simulation study.

 

 

Please, also, how do I write it as vector?

For example:

seq( X__B(X__0, SX__1[k], SX__2[k]), k=1..N);

HFloat(2.2461753368121835), HFloat(2.386372343185135), HFloat(1.305439309043902), HFloat(1.3279333023039865), HFloat(1.2430844367475704)

 

I copied the comand and the results from your post but the results are not listed in vector form; I mean not contain block brackets [ ] 

 

Thank you.

 

SECOND SCENARIO...

Josolumoh 25
March 16 2019
0 0

@tomleslie 

Thank you Prof,

 

The second scenario discribes by interests:

1. I want to sample random numbers for X_1 and X_2. Afterwards, I will export the sampled vectors and keep

2. I want to use the vectors to evaluate X_B (Yes! " the function X_B does not define a new random variable") but results from substituting X_0, X_1, and X_2. 

3. These I think the command attached will sove. However, my attempt to run the copied command produced  " Error, invalid arrow procedure"

 4. Please, help me correct it.

 

My regards.


 

> 

restart; with(Statistics); randomize(); B := unapply*(simplify*binomial(x+r-1, x)*((1/2)*b/(d*x+1+(1/2)*b))^r*((d*x+1)/(d*x+1+(1/2)*b))^x/(d*x+1), [r, b, d]); beta := Array(0 .. 2, [.25, .6, -.2]); N := 5

155277071509912

  

unapply*(simplify*binomial(x+r-1, x)*((1/2)*b/(d*x+1+(1/2)*b))^r*((d*x+1)/(d*x+1+(1/2)*b))^x/(d*x+1), [r, b, d])

  

beta := Array(0..2, {(1) = .25, (2) = .6})

  

5

 (1)
> 

X__0 := 1; X__1 := RandomVariable(Bernoulli(1/2)); X__2 := RandomVariable(Normal(0, 1)); SX__1 := Sample(X__1, N); SX__2 := Sample(X__2, N)

`#msub(mi("X"),mi("0"))` := 1

  

`#msub(mi("X"),mi("1"))` := _R11

  

`#msub(mi("X"),mi("2"))` := _R12

  

`#msub(mi("SX"),mi("1"))` := Vector[row](5, {(1) = 1.0, (2) = .0, (3) = .0, (4) = 1.0, (5) = 1.0}, datatype = float[8])

  

`#msub(mi("SX"),mi("2"))` := Vector[row](5, {(1) = -0.9622324619413775e-1, (2) = 1.9771608385550665, (3) = 1.6198494896183102, (4) = -.49236514936140713, (5) = -0.25599340478732677e-1}, datatype = float[8])

 (2)
> 

"#`X__B:= `(x, y, z)-> exp(add(beta[k]*[x,y,z][k+1], k = 0 .. 2)); `X__B`:= (x, y, z)-> local k; exp ( add ( beta[k]*[x,y,z][k+1], k = 0 .. 2 ) ); "

Error, invalid arrow procedure

"#`X__B:= `(x, y, z)-> exp(add(beta[k]*[x,y,z][k+1], k = 0 .. 2)); X__B:= (x, y, z)-> local k; exp ( add ( beta[k]*[x,y,z][k+1], k = 0 .. 2 ) ); "

  
> 

``


 

Download Second_Scenario.mw

 

Re: Wrong /Incorrect Estimates...

Josolumoh 25
March 16 2019
0 0

@mmcdara 

 

Please, help me out on that problem.

 

Wrong Estimates...

Josolumoh 25
March 15 2019
0 0

@mmcdara 

 

Hello Please,

 

I discovered that when subtituting the values of X_0,  X_1 and X_2 generated to exp(0.25+0.6*_R -0.2_R0)

the results are different to what the computer generated.

For example:

From the results below;

X_0 = 1

X_1 =Vector[row]([1., 1., 0., 0., 0.])

X_2 =Vector[row]([.203855402371112, -0.988718240271234e-1, -0.826980970419822e-1, -.168119128121408, .162021299867933])

Consequently,

exp(.25+.6*_R-.2*_R0)

Sample_b := Vector[row](5, {(1) = .9726008729365625, (2) = 2.150871828664164, (3) = 1.4046239626237194, (4) = 1.4658659638712146, (5) = 2.07645686621087}, datatype = float[8])

But, when subtituted manually, I get different results as we can see below

1.  exp(.25+.6*1-.2*0.203855402371112)  =  2.246175
2,  exp(.25+.6*1-.2*-0.988718240271234e-1) = 2.386372

3.  exp(.25+.6*0-.2*-0.826980970419822e-1) = 1.305439

4.  exp(.25+.6*0-.2*-.168119128121408) =  1.327933

5. exp(.25+.6*0-.2*0.162021299867933) = 1.243084

 

Please, how do I make (exp(.25+.6*_R-.2*_R0)) use the values of  X_0,  X_1, and X_2 generated and return exact answers instead of it generating another ones?

> 

restart
> 

> 

with(Statistics):
> 

B := unapply(simplify(binomial(x+r-1, x)*((1/2)*b/(1+d*x+(1/2)*b))^r*((d*x+1)/(1+d*x+(1/2)*b))^x/(d*x+1)), [r, b, d]);

proc (r, b, d) options operator, arrow; binomial(x+r-1, x)*(b/(2*d*x+b+2))^r*2^x*((d*x+1)/(2*d*x+b+2))^x/(d*x+1) end proc

 (1)
> 

beta := Array(0 .. 2, [.25, .6, -.2]);

beta := Array(0..2, {(1) = .25, (2) = .6})

 (2)
> 

N := 5:
> 

X__0 := 1:
> 

X__1 := RandomVariable(Bernoulli(1/2));

_R

 (3)
> 

Sample(X__1, N)

Vector[row]([1., 1., 0., 0., 0.])

 (4)
> 

X__2 := Statistics:-RandomVariable(Normal(0, 1));

_R0

 (5)
> 

Sample(X__2, N)

Vector[row]([.203855402371112, -0.988718240271234e-1, -0.826980970419822e-1, -.168119128121408, .162021299867933])

 (6)
> 

X__B := exp(add(`~`[`*`](beta[k], X__ || k), k = 0 .. 2));

exp(.25+.6*_R-.2*_R0)

 (7)
> 

Sample_b := Sample(X__B, N)

Sample_b := Vector[row](5, {(1) = .9726008729365625, (2) = 2.150871828664164, (3) = 1.4046239626237194, (4) = 1.4658659638712146, (5) = 2.07645686621087}, datatype = float[8])

 (8)
> 

``


 

Download Wrong_Estimates_of_X_B.mw
 

> 

restart
> 

> 

with(Statistics):
> 

B := unapply(simplify(binomial(x+r-1, x)*((1/2)*b/(1+d*x+(1/2)*b))^r*((d*x+1)/(1+d*x+(1/2)*b))^x/(d*x+1)), [r, b, d]);

proc (r, b, d) options operator, arrow; binomial(x+r-1, x)*(b/(2*d*x+b+2))^r*2^x*((d*x+1)/(2*d*x+b+2))^x/(d*x+1) end proc

 (1)
> 

beta := Array(0 .. 2, [.25, .6, -.2]);

beta := Array(0..2, {(1) = .25, (2) = .6})

 (2)
> 

N := 5:
> 

X__0 := 1:
> 

X__1 := RandomVariable(Bernoulli(1/2));

_R

 (3)
> 

Sample(X__1, N)

Vector[row]([1., 1., 0., 0., 0.])

 (4)
> 

X__2 := Statistics:-RandomVariable(Normal(0, 1));

_R0

 (5)
> 

Sample(X__2, N)

Vector[row]([.203855402371112, -0.988718240271234e-1, -0.826980970419822e-1, -.168119128121408, .162021299867933])

 (6)
> 

X__B := exp(add(`~`[`*`](beta[k], X__ || k), k = 0 .. 2));

exp(.25+.6*_R-.2*_R0)

 (7)
> 

Sample_b := Sample(X__B, N)

Sample_b := Vector[row](5, {(1) = .9726008729365625, (2) = 2.150871828664164, (3) = 1.4046239626237194, (4) = 1.4658659638712146, (5) = 2.07645686621087}, datatype = float[8])

 (8)
> 

``


 

Download Wrong_Estimates_of_X_B.mw

 

Code posted...

Josolumoh 25
March 05 2019
0 0


 

 

Hello please,

 

Here is the code.
The code is taking very long time to run and some times it stops working, maybe because of large numbers involve.

Is there a way of making it work faster? NULL

 

My regards 


 

Download Random_Number_Generation_from_QNBR3.mw
 

restart*with(Statistics); B := unapply(simplify(binomial(x+r-1, x)*((1/2)*b/(1+d*x+(1/2)*b))^r*((d*x+1)/(1+d*x+(1/2)*b))^x/(d*x+1)), [r, b, d]); beta := Array(0 .. 2, [.25, .6, -.2]); N := 1000; X__0 := 1; X__1 := RandomVariable(Bernoulli(1/2)); X__2 := RandomVariable(Normal(0, 1)); X__B := exp(add(`~`[`*`](beta[k], X__ || k), k = 0 .. 2)); Sample_b := Sample(X__B, N); SUBS := proc (n) options operator, arrow; [r = 5, b = Sample_b[n], d = .2] end proc; Xmax := 1000; NperB := 10000; U := RandomVariable(Uniform(0, 1)); AllSamples := NULL; for n to N do PointwiseB := [seq([x, B(`~`[rhs](SUBS(n))[])], `in`(x, [`$`(0 .. Xmax)]))]; ShouldBe1 := add(`~`[op](2, PointwiseB)); PointwiseB := [seq([x, B(`~`[rhs](SUBS(n))[])/ShouldBe1], `in`(x, [`$`(0 .. Xmax)]))]; c := CumulativeSum(`~`[op](2, PointwiseB)); ApproximateCDF := zip(proc (u, v) options operator, arrow; [u[1], v] end proc, PointwiseB, convert(c, list)); AllSamples := AllSamples, map(proc (u) options operator, arrow; ListTools:-BinaryPlace(`~`[op](2, ApproximateCDF), u) end proc, Sample(U, NperB)) end do; AllSamples := `<,>`(AllSamples); ExportMatrix("C:/Users/jamiu.olumoh/Documents/QNBRsimfiles/Sample1000502.csv", AllSamples, target = csv, format = rectangular, mode = ascii); M[1] := Matrix(1000, 1, datatype = float[8]); S[1] := Statistics:-Sample(X__1, M[1]); ExportMatrix("C:/Users/jamiu.olumoh/Documents/QNBRBeta1/X11.csv", S[1], target = csv); M[2] := Matrix(1000, 1, datatype = float[8]); S[2] := Statistics:-Sample(X__2, M[2]); ExportMatrix("C:/Users/jamiu.olumoh/Documents/QNBRBeta2/X21.csv", S[2], target = csv); for n to N do PointwiseB := [seq([x, B(`~`[rhs](SUBS(n))[])], `in`(x, [`$`(0 .. Xmax)]))]; ShouldBe1 := add(`~`[op](2, PointwiseB)); pl || n := plot(B(`~`[rhs](SUBS(n))[])/ShouldBe1, x = 0 .. 5) end do; plots:-display(seq(pl || n, n = 1 .. N))

  

NULL


 

Download Random_Number_Generation_from_QNBR3.mw

 

@tomleslie 

Re: Program running slow...

Josolumoh 25
March 04 2019
0 0

@tomleslie 

 

Please help me check the program, if there is anything wrong with it.

Or is there  anything to be added to make it work faster ?

Thank you please.

Program running slow...

Josolumoh 25
March 01 2019
0 0

@mmcdara

Please help me check the programe as it is taken very long time to complete a process and at time stops wlthout completing the task.

Anything to be added to make run fast?

 

Thank You...

Josolumoh 25
February 27 2019
0 0

@tomleslie 

 

 

 

Accept my appreciation please.

Export to csv...

Josolumoh 25
February 27 2019
0 0

@tomleslie 

Wow!

 

This idea appears to be easier but still not working. Please check it. 

 

resarts;

N: =50:

X_1 := RandomVariable(Bernoulli(1/2));

X_2 := RandomVariable(Normal(0,1));

ExportVector( "C:/Users/jamiu.olumoh/Desktop/testvecX1.csv", X_1, target=csv);

ExportVector( "C:/Users/jamiu.olumoh/Desktop/testvecX2.csv", X_2, target=csv);

 

 

 

 

 

 

Thank You...

Josolumoh 25
February 26 2019
0 0

@acer 

 

 

Many thanks. It works...

Exporting as csv...

Josolumoh 25
February 26 2019
0 0

@tomleslie 

 

Now you understand my question.

I wish to export it directly as csv. Like exporting the matrix, the command below works well but not for the vector

    "ExportMatrix("C:/Users/jami.olu/Documents/X1.csv", X__1, target = csv)"  

 

 

1 2 3 4 Page 1 of 4