Maple programs

Rolling of a die

This program simulates the rolling of a die times.

> DP:=proc(n) local i,v,die;

DP is the name of the program, are the brackets inside which the program is located. At the beginning of the program we listed all the variables used in it.

> die:=rand(1..6);

That defines as a randomly chosen integer between 1 and 6.

> v:=[0,0,0,0,0,0];

We'll use array for counting of the appearances of 1, 2, ..., 6.

> to n do i:=die(); v[i]:=v[i]+1 od;

We roll a die times and every time add 1 to the corresponding counter.

> print(v) end:

We print the final numbers of the occurrences of 1, 2, ..., 6 after rolling a die times.

> to 5 do DP(6000) od;

We executed our program 5 times with every time.

Buffon's needle

I was told that the program BuffonsNeedle accompanying the book, doesn't work. To correct it, add one line:

:

This program does the same thing as the BuffonsNeedle, only without plotting.

> BN:=proc(n) local d,a,i; i:=0;

BN is the name of the program, is the number of times we drop a needle, are the brackets inside which the program is located. At the beginning of the program we listed all the variables used in it and initialized variable counting the number of intersections.

> to n do d:=evalf(rand()/10^12); a:=evalf(rand()/10^12*Pi/2); if(d<sin(a)) then i:=i+1 fi od;

We drop a needle times adding 1 to if it intersects with a line.

> print(evalf((n+n)/i)); end:

We print the approximation to given by this simulation.

> to 5 do BN(1000) od;

We executed our program 5 times with every time.

Monte Carlo method

The program MonteCarlo accompanying the text, evaluates the probability that a randomly chosen point is lying under the curve instead of the evaluating of the area under the curve. It is easy to correct that changing print(evalf(counter/n)); at the end of the program to print(evalf(counter/n*(xmax-xmin)*ymax)); .

This program evaluates the area between two curves, and from to using the Monte Carlo method with points. Functions must be entered with arrows, like x->x^3, for instance.

> MC:=proc(f,g,a,b,n) local c,i,MIN,MAX,p,pf,pg,r,v,vf,vg;c:=0;

MC is the name of the program, are the brackets inside which the program is located. At the beginning of the program we listed all the variables used in it and initialized the counter .

MIN:=min(f(a),f(b),g(a),g(b)); MAX:=max(f(a),f(b),g(a),g(b));

We choose the minimal and maximal values of our functions at the ends of the segment as initial approximations to the absolute minimum and maximum.

> for i to n do v[i]:=evalf(a+rand()/10^12*(b-a));

We randomly choose points on the interval from to .

> vf[i]:=evalf(f(v[i])); vg[i]:=evalf(g(v[i]));

We evaluated our functions in these randomly chosen points.

> MIN:=min(MIN,vf[i],vg[i]); MAX:=max(MAX,vf[i],vg[i]) od;

If some of the values are less then the previous minimum, or more than the previous maximum, we decrease minimum, or increase maximum correspondingly.

> for i to n do r[i]:=evalf(MIN+rand()/10^12*(MAX-MIN));

Now we randomly choose the -coordinates between minimal and maximal values of our functions for all previously chosen random values of -coordinate.

> if ((vf[i]<=r[i] and r[i]<vg[i]) or (vg[i]<=r[i] and r[i]<vf[i])) then c:=c+1 fi od;

If the point is lying between graphs of our functions, we add 1 to the counter .

> print(evalf((b-a)*(MAX-MIN)*c/n));

We print the estimated area.

> with(plots):

This line starts the package . We need it to use and below.

> p:=pointplot({seq([v[i],r[i]],i=1..n)}):

We prepare for drawing our points.

> pf:=plot(f(x),x=a..b,y=MIN..MAX):

This is the graph of .

> pg:=plot(g(x),x=a..b,y=MIN..MAX):

This is the graph of .

> display(p,pf,pg)

This command displays the graphs of the functions and the points on one plot.

> end:

The end of the program.

> for k to 5 do print(n=200*k); MC(x->x,x->x^3,-1,1,200*k) od;

We run our program 5 times for functions and on the interval from -1 to 1 using 200, 400, 600, 800 and 1000 points.

> evalf(int(abs(x^3-x),x=-1..1));

This is the actual area between the curves.

Binomial and Poisson Distributions

> BP:=proc() local b, n, p, k, Po;

interface(prompt="Enter the number of trials: "); n:=parse(readline(),statement); interface(prompt="Enter the probability of success: "); p:=parse(readline(),statement); interface(prompt="Enter the number of successes: "); k:=parse(readline(),statement);

interface(prompt="> "); b:=(n,p,k)->binomial(n,k)*p^k*(1-p)^(n-k); Po:=(lambda,k)->lambda^k/k!*exp(-lambda); print(Binomial=evalf(b(n,p,k)), Poisson=evalf(Po(n*p,k)))

end:

> BP();

Enter the number of trials: 1000

Enter the probability of success: 0.001

Enter the number of successes: 5

Leading Digits

> LD:=proc(f,n) local i,v,L,w; L:=k->floor(f(k)/10^ilog10(f(k))); v:=[0,0,0,0,0,0,0,0,0]; w:=v; for i to n do v[L(i)]:=v[L(i)]+1 od; print(Results=v); for i to 9 do w[i]:=round(log[10](1.+1/i)*n) od; print(Benford=w) end:

> LD(k->2^k,100);

> LD(k->3*2^k,100);

> LD(k->2^k,1000);

> LD(k->3*2^k+2*3^k,500);

> with(combinat):

> LD(k->fibonacci(k),1000);

Second Digits

> SD:=proc(f,n) local i,j,g,L,v,w,m; L:=(f,k)->floor(f(k)/10^ilog10(f(k))); g:=k->f(k)-L(f,k)*10^ilog10(f(k));v:=[0,0,0,0,0,0,0,0,0,0]; w:=v; for i to n do j:=`if`(ilog10(g(i))=ilog10(f(i))-1,L(g,i)+1,1); v[j]:=v[j]+1 od; print(Results2=v); for i to 10 do w[i]:=round(n*sum(log[10](1.+1/(10*m+i-1)),m=1..9)) od; print(Benford2=w) end:

> SD(k->2^k,100);

> SD(k->3*2^k,100);

> SD(k->2^k,1000);

> SD(k->3*2^k+2*3^k,500);

> with(combinat):

> SD(k->fibonacci(k),1000);

N Fold Convolution

> NFC:=proc(g,n) local c, cred, i, m, p, t, v, Points, NF; p:=g^n; t:=indets(g)[1]; c:=seq([i,coeff(p,t,i)],i=ldegree(p)..degree(p)); cred:={c}; for i to degree(p)-ldegree(p)+1 do if c[i,2]=0 then cred:=cred minus {c[i]} fi od; with(plots); Points:=pointplot(cred); m:=subs(t=1,diff(g,t)); v:=subs(t=1,diff(g,t$2))+m-m^2;NF:=plot(exp(-(x-m*n)^2/(2*n*v))/sqrt(2*Pi*n*v),x=ldegree(p)..degree(p)); display(Points,NF) end:

> NFC(.25*t+.25*t^2+.5*t^3,10);

> NFC(.25*t+.25*t^2+.5*t^3,200);

> NFC(.25*t+.75*t^3,10);

> NFC(.25*t+.75*t^3,100);

> NFC(.25*t+.75*t^3,200);

> NFC(.25*t+.01*t^2+.74*t^3,10);

> NFC(.25*t+.01*t^2+.74*t^3,100);

> NFC(.25*t+.01*t^2+.74*t^3,200);

95% Confidence Interval for Bernoulli Trials

This is the solution of Excercise 9.1.11 (p. 339).

> CB:=proc(prob,n) local N,p,q,r; N:=0; to n do r:=rand()/10^12; if r<prob then N:=N+1 fi od; p:=N/n; q:=1-p; evalf([p-2*sqrt(p*q/n),p+2*sqrt(p*q/n)]) end:

> CB(.3,10000);

> k:=0: to 100 do v:=CB(.3,10000): if v[1]<0.3 and 0.3<v[2] then k:=k+1 fi od: k;