 # ranxi() gives a random # uniform on (0,1), using generic maple RNG  
ranxi := () -> (10e-13)*rand();            

 # ranexp(lambda) should give a random # distributed according to the 
 #    exponential probability density function exp[-x/lambda]/lambda
 #    where x is on ( 0, +infinity ) and lambda ( >0 ) is a parameter.   
ranexp := (lambda) -> ranxi()            ;
 #                    ^^^^^^^^^^^^^^^^^^^ <-- MODIFY -- fill in correct form
####
 # program that calls ranexp(lambda) to generate a histogram of 
 #    distribution and also measures the first moment <x>
run := proc(lambda, n_mea, n_smpls_per_mea)

global xacc, hist, bin_range, bin_size, n_bins, bin_size_inv, midpt, bin;
local bi, st, m, i, x;
 # -- define parameters and procedures for histogram:
bin_range := 6. *lambda;                #NOTE: prefactor can be changed
bin_size := 0.1 *lambda;                #NOTE: prefactor can be changed
n_bins := trunc(bin_range/bin_size)+1;  #last bin to store all x > bin_range
bin_size_inv := 1./bin_size;
hist := array(1..n_bins);
for bi to n_bins do hist[bi] := 0 od:
bin  := proc(x)    
    global n_bins, bin_size_inv; 
    min(trunc(x*bin_size_inv)+1,n_bins) 
end:                                    # bin(x) gives which bin x is in.  
midpt := proc(i) global bin_size; (i-0.5)*bin_size end:
                                        # midpt(i) gives mid-pt of bin i

 # -- xacc[m] stores the m^th measurement of <x>
xacc := array(1..n_mea);

 # -- n_mea measurements [first do loop (with m)], 
 #    EACH with n_smpls_per_mea samples [second do loop (with i)] :
for m from 1 to n_mea do
    xacc[m] := 0:
    for i to n_smpls_per_mea do
        x := ranexp(lambda):
        xacc[m] := xacc[m] + x:
        bi := bin(x):
        hist[bi] := hist[bi] +1
    od;
    xacc[m] := xacc[m]/n_smpls_per_mea:
od;

end:
