## EXERCISE 1
_Christian Petersen (s164568)_
First we experimented with different values for **a** **c** and **m**.

We found that the best values where when M was an exponential of 2 and **c** > **a** while **a** should always be odd.

I.e


So when we choose **a = 4** (even) then we see a very skewed value.

When we do a scatter plot of $U_i$ vs $U_{i+1}$ for our generated LCG distribution:


It does look like is a lot of correlation here, first of all the we would expect there to be $C(16,2)=105$ dots if it was completely uncorrelated, i.e all the possible pairs would be chosen but this is obviously not the case, we only have 8 dots on this plot.

### $\chi^2$ Test

def chi_squared(hist):
  	#hist[0] is the count in each bin/class, it is of length 10.
    leng = len(hist[0])
    # we expect, in a uniform distribution the amount of datapoints distributed evenly:
    expected = np.sum(hist[0])/leng
    # create the list for T
    l = [((n - expected)**2)/expected for n in hist[0]]

    chi = sum(l)
    return chi

We get a $\chi = 3.686$ which corresponds to a $p = 0.9026$ with a $df = 9$ for 10 bins

### K-S Test

For the K-S test we get the following values
$D_n = 0.93$ So what we cannot disregard the null hypothesis. In other words the sample data does indeed follow that of a uniform distribution.
### Runtest
#### Runtest 1 Conradsen (Above Below)

$H_0$ : the sequence was produced in a random matter
$H_1$ : the sequence is not random

We compute with $\mathcal{N}\left(2\cdot \dfrac{n_1\cdot n_2}{n_1+n_2}+1, 2\dfrac{n_1 n_2(2n_1 n_2-n_1 -n_2)}{(n_1+n_2)^2(n_1 + n_2-1)}\right)$

$=\mathcal{N}(5001,2499.7)$

Then we find the Z value.
$Z = \frac{n-\mu}{\sigma}$
where is the total number of runs.
$Z = \frac{3125-5001}{\sqrt{2499.7}}=-37.52$
Since our z value is not between -1.96 and 1.96 we cannot reject $H_0$

## EXERCISE 2
_Christian Petersen (s164568)_
To generate our Geometric distribution I used the numpy library within python:

xs = np.random.uniform(0, 1)

Then I generated the geometric sequence the following way:

p = 0.1
size = 100000
i = 0
binom_dist = []
while len(binom_dist) <= size:
    i +=1
    x = np.random.uniform(0,1)
    if x > p:
        print("fail")
    else:
        print("success")
        binom_dist.append(i)
        i = 0


For the test we chose our _p_ value as _p = 0.1_ and generated 10 000 outcomes. The mean of the simulation above is 9.9173, so the _emperical p_ value is = 0.092. Overall this is a pretty close outcome and that is expected with such a high sampling size.

If we combine this with the theoretical values to see how close they are to each other we see that the emperical distribution does indeed follow the theoretical.


Now to simulate the following 6 point distribution:
We will generate the distribution using the 3 methods.

1. Direct (crude) method
2. Rejection method
3. Alias method

### Direct method

def crude_method(distribution):
    k = len(distribution)
    cdf = np.cumsum(distribution)
    U = np.random.uniform()

    for i in range(k):
        val = cdf[i]
        if U <= val:
            return i+1
        else:
            print("U : ", U)
            print("next, ", val)
    return i+1

#### Results
After 1000 runs

### Rejection method

def rejection_method(distribution):
    k = len(distribution)
    cdf = np.cumsum(distribution)
    c = 0.8

    while True:
        U1 = np.random.uniform()
        U2 = np.random.uniform()
        I = 1 + k*U1
        I_int = int(I)-1
        print(I_int)

        if U2 < distribution[I_int]/c:
            return I_int

Results for 1000 runs:
_Note the graph starts from 0 instead of 1_

### Alias Method
#### Results

## EXERCISE 3
_Nicholas Rose (s164580)_

We begin with generating simulated values from the exponential distribution.

def exp(lam, n):
    data = []
    while len(data) < n:
        u = random.uniform(0,1)
        data.append(-math.log(u)/lam)
        
    return data

The function takes arguments lambda (the mean) and n (the desired number of simulated values). With lambda = 1.5 and n = 10000 (a sufficiently large number to generate a reliable trend), the histogram displays the following trend.


The next simulated distribution will be the Normal Distribution with Box-Muller transformation.

def box(n):
    data = []
    while len(data) < n:
        u1 = random.uniform(0,1)
        u2 = random.uniform(0,1)
        z1 = math.sqrt(-2*math.log(u1)) * math.cos(2 * math.pi * u2)
        z2 = math.sqrt(-2*math.log(u1)) * math.sin(2 * math.pi * u2)
        data.append(z1)
        data.append(z2)
    return data

We assume a mean of 0 and generate n = 10000 simulated data points. The histogram resulting from this simulation is displayed below.


The final distribution that we simulated is the Pareto distribution with $\beta$ = 1 and numerous different k values.

def paret(n,k):
    data = []
    while len(data) < n:
        u = random.uniform(0,1)
        data.append((u**(-1/k))-1)
        
    #[beta, infinity]
    mean = k/(k-1)
    var = k/(((k-1)**2) * (k-2))
    
    return data, mean, var
             

With k = 2.5, the resulting histogram appears as so:


mean = 1.6666666666666667
var = 2.2222222222222223

With k = 4:


1.3333333333333333,
0.2222222222222222)

With k = 10:

mean = 1.1111111111111112
var = 0.015432098765432098

As we increase the k, the distribution variance decreases and our simulated values become more representative of the real Pareto distribution.

##CONTINUED
Using the Box-Muller method for the Normal distribution, we generate 10 data sets of 10 variables and use this to generate 10 different values of the mean, variance and confidence intervals displayed below.

| mean | variance | lower | upper |
|------------- |------------ |------------- |------------- |
| -0.08613556 | 0.15556972 | -0.978383 | 0.80611188 |
| 0.01617023 | 0.05422704 | -0.51061154 | 0.542952 |
| 0.02171605 | 0.14317181 | -0.83424015 | 0.87767224 |
| 0.02955248 | 0.20349708 | -0.99092132 | 1.05002629 |
| -0.41329691 | 0.14034127 | -1.26074963 | 0.43415581 |
| 0.25354933 | 0.07046481 | -0.34694502 | 0.85404367 |
| -0.34476772 | 0.10396389 | -1.07416483 | 0.38462939 |
| -0.04947222 | 0.107471 | -0.79106998 | 0.69212555 |
| -0.01323924 | 0.03520109 | -0.43766413 | 0.41118566 |
| 0.03099813 | 0.0452292 | -0.45009841 | 0.51209468 |

## EXERCISE 4
_Christian Petersen (s164568)_

when _n = 10_ and _mst_ = 8, _mct_ =1

We can service 0.8 per unit.

Initially we use an exponential distribution with a mean of 1 to model the service time. In the simulation code `arrival_model()`and `service_time()`are functions that sample from a exponential distributions

The simulation code is defined in the following main function:

def main(offered_traffic, units):
    buffer = []
    blocked_count = 0
    # put one customer in buffer
    counter = 0
    arrival_t = arrival_model()
    service_t = service_time()
    buffer.append(service_t)
    # init
    while offered_traffic > counter:
        if len(buffer) > 0:
            index_min = min(range(len(buffer)), key=buffer.__getitem__)
        else:
            buffer.append(service_time())
            index_min = min(range(len(buffer)), key=buffer.__getitem__)

        buffer.sort()
        # print("Pre:  ", buffer)
        # print("  ")
        # arrival handler

        # if arrival time is less than minimum service time
        if arrival_t < buffer[index_min]:
            # print("   ")
            # print("arrival handler ",)
            # print("")
            # print("arrival_t = ", arrival_t, "service_t = ", buffer[index_min])

            # if the buffer is full:
            if len(buffer) >= units:
                # print("buffer of size: ",len(buffer), "greater than ",units)
                # add one to block counter
                blocked_count += 1
                # subtract arrival time from buffer
                for i in range(len(buffer)):
                    buffer[i] = buffer[i] - buffer[index_min]
                buffer.pop(index_min)

            else:
                # designate service time
                service_t = service_time()
                buffer.append(service_t)
                buffer.sort()
                # generate new arrival timew
                for i in range(len(buffer)-1):
                    #print("doing subtraction", buffer[i],"-", arrival_t)
                    buffer[i] = buffer[i] - arrival_t
                arrival_t = arrival_model()

        # service handler
        elif buffer[index_min] < arrival_t:
            # print("")
            # print("service handler: ")
            # print("")
            for i in range(len(buffer)):
                buffer[i] = buffer[i] - buffer[index_min]
            buffer.pop(0)
        counter +=1
        # print("Post: ",buffer)
    return blocked_count/offered_traffic

In order to calculate the confidence interval for our blocked customer ratios (BCRs) we run our program 10 times with offered traffic being 10 000

$BCRs= [0.1319, 0.1188, 0.1228, 0.1264, 0.1286, 0.1266, 0.1465, 0.1267, 0.1334, 0.1316]$


And a 95% confidence interval:

`[0.12970850969659872, 0.13955149030340133]`

We repeat the process for a gamma distribution:

$BCRs = [0.1367, 0.1238, 0.1267, 0.1393, 0.1311, 0.107, 0.1245, 0.1209, 0.1265, 0.0967]$

$CI = [0.11408818124127314, 0.13255181875872687]$

And for the hyperexponential distribution:

$BCRs = [0.1526, 0.1446, 0.1638, 0.1481, 0.1389, 0.1411, 0.1532, 0.1452, 0.1289, 0.1558]$

$CI =[0.1401978489320671, 0.15424215106793288]$


Now for experimenting with the service time, we tried a constant service time of 8 units.


We can now see that the blocking rate is much lower now.

Finally we tried with a normal distribution with mean 8 and variance 1.

$BCRs = [0.1534, 0.1571, 0.1701, 0.1501, 0.1555, 0.1639, 0.1539, 0.1682, 0.163, 0.1]$

$CI = [0.13925067474568234, 0.16778932525431764]$

## EXERCISE 5

_Nicholas & Jack Rose (s164580,s164559)_

In this exercise we estimate the result for the integral of $e^x$ between 1 and 0. The actual value of this is 1.71828182846.
### Monte Carlo Estimator
For the Monte Carlo estimator we first generate an array of 100 random numbers between 0 and 1 with a uniform distribution. We then sum the the result of the exponential function to the power of these random numbers and find the average. This gives us our final estimate.

def montcarlo(samples):
    randno = np.random.uniform(0, 1, samples)
    expo = lambda t: math.exp(t)
    func = np.vectorize(expo)
    sum = 0
    for i in range(samples):
        sum += math.exp(randno[i])
    return sum/samples, func(randno)```

$mean = 1.7877237743041818$
$95 % Confidence interval :   low bound = 1.6334373079992377   high bound = 1.8257389573189353

As you can see this result is reasonably close to the actual value.

### Antithetic variables
Antithetic variables make the estimate of our result a lot more accurate and it is also a lot more efficient because each random variable is used twice. 

def antithetic(samples):
randno = np.random.uniform(0, 1, samples)
expo = lambda t: (math.exp(t) + (math.exp(1)) / math.exp(t))/2
func = np.vectorize(expo)
sum = 0
for i in range(samples):
sum += (math.exp(randno[i]) + (math.exp(1)) / math.exp(randno[i]))/2
return sum/samples, func(randno)

$mean = 1.7268840302242356$
$95 % Confidence interval :   low bound = 1.714354100428142  high bound = 1.739413960020328

Here it shows that by using this method we have a more narrow confidence interval and a more accurate average.

### Control Variates

def controlvariate(samples):
randnom = np.random.uniform(0, 1, samples)
expo = lambda t:math.exp(t)
func = np.vectorize(expo)
randnot = func(randnom)
summ = 0
sumt = 0
summt = 0
meanm = randnom.mean()
meant = randnot.mean()
for i in range(samples):
summ += (randnom[i] - meanm)**2
sumt += (randnot[i] - meant)**2
summt +=(randnom[i] - meanm)*(randnot[i] - meant)
varm = summ/(samples-1)
vart = sumt/(samples-1)
cov =summt/(samples-1)
return vart - (cov**2)/varm, randnot

$Variance = 0.0038754777181509437   

### Stratified Sampling

def stratifiedsampling(samples,sampleselection):
randno = np.random.uniform(0, 1, samples)
rand1 = []
rand2 = []
rand3 = []
rand4 = []
rand5 = []
rand6 = []
rand7 = []
rand8 = []
rand9 = []
rand10 = []
for i in range(samples):
if randno[i] <= 0.1:
rand1.append(randno[i])
elif randno[i] <= 0.2:
rand2.append(randno[i])
elif randno[i] <= 0.3:
rand3.append(randno[i])
elif randno[i] <= 0.4:
rand4.append(randno[i])
elif randno[i] <= 0.5:
rand5.append(randno[i])
elif randno[i] <= 0.6:
rand6.append(randno[i])
elif randno[i] <= 0.7:
rand7.append(randno[i])
elif randno[i] <= 0.8:
rand8.append(randno[i])
elif randno[i] <= 0.9:
rand9.append(randno[i])
else:
rand10.append(randno[i])

wi = (math.exp(rand1[sampleselection]/10) + math.exp(1/10 + rand2[sampleselection]/10) + math.exp(2/10 + rand3[sampleselection]/10)
+ math.exp(3/10 + rand4[sampleselection]/10) + math.exp(4/10 + rand5[sampleselection]/10) + math.exp(5/10 + rand6[sampleselection]/10)
+ math.exp(6/10 + rand7[sampleselection]/10) + math.exp(7/10 + rand8[sampleselection]/10) + math.exp(8/10 + rand9[sampleselection]/10)
+ math.exp(9/10 + rand10[sampleselection]/10))/10
return wi, randno

$mean = 1.7326506167199665$

We use stratified sampling when we already have some knowledge of the data and are able to split it up into distinct categories. From this we then select a sample at a chosen index and find the estimate from these chosen random numbers. As we know that the data has a small value for the covariate there will be a smaller error in the estimation.


## EXERCISE 6
_Jack Rose (s164558)_

In the random walk Metropolis Hastings algorithm we only accept a state, which is generated by randomly adding or taking one from the previous state, if the value of g(x)<=g(y) where x is the current state and y is the next state. The state is also accepted if g(y)/g(x) is greater than a random uniform number between 0 and 1. This allows the data to include some of the outliers of the distribution making it a more accurate representation of the actual distribution of the formula. 

d
ef met_hast_alg(startval):
list_ = [startval]
for i in range(10000):
x= list_[i]
y = x + random.choice([-1, 1])
if y < 0:
y = 10
if y > 10:
y = 0
g1 = (8**x)/math.factorial(x)
g2 = (8**y)/math.factorial(y)
if g2 >= g1:
list_.append(y)
elif g2<g1:
u = np.random.rand()
if u < g2/g1:
list_.append(y)
else:
list_.append(x)
return list_

Here we can see the predicted distribution of the formula.
![](https://i.imgur.com/IUwPBle.png)

Here we can see the actual distribution of the formula
![](https://i.imgur.com/MAvesLW.png)

For the coordinate and direct Metropolis Hastings algorithm we simulate the distribution of a function with two variables. The boundary for these variables is that 0 is less than or equal to i + j which is less than or equal to 10. were implemented as shown below

def met_hast_alg_coowise(startval1,startval2):
list_ = [np.array([startval1, startval2])]
for i in range(10000):
x1 = list_[i][0]
y1 = list_[i][1]
rand1 = random.choice([-1, 1])
if rand1 == -1:
y2 = y1
x2 = x1 + random.choice([-1, 1])
if x2 + y2 > 10:
x2 = 0
if x2 < 0:
x2 = 10 - y2
elif rand1 == 1:
x2 = x1
y2 = y1 + random.choice([-1, 1])
if x2 + y2 > 10:
y2 = 0
if y2 < 0:
y2 = 10 - x2
g1 = (4 ** x1 / math.factorial(x1)) * (4 ** y1 / math.factorial(y1))
g2 = (4 ** x2 / math.factorial(x2)) * (4 ** y2 / math.factorial(y2))
if g2 >= g1:
list_.append(np.array([x2,y2]))
elif g2 < g1:
u = np.random.rand()
if u < (g2 / g1):
list_.append(np.array([x2,y2]))
else:
list_.append(np.array([x1,y1]))
return list_

def met_hast_alg_dir(startval1,startval2):
list_ = [np.array([startval1, startval2])]
for i in range(10000):
x1 = list_[i][0]
y1 = list_[i][1]
rand1 = random.choice([-1,0,1])
if rand1 == -1:
y2 = y1
x2 = x1 + random.choice([-1, 1])
if x2 + y2 > 10:
x2 = 0
if x2 < 0:
x2 = 10 - y2
elif rand1 == 0:
rand2 = random.choice([-1, 1])
rand3 = random.choice([-1, 1])
x2 = x1 + rand2
y2 = y1 + rand3
if x2 + y2 > 10:
if rand2 == 1 & rand3 < 1:
x2 = 0
elif rand2 < 1 & rand3 == 1:
y2 = 0
elif rand2 == 1 & rand3 == 1:
y2 = 0
if x2 + y2 <0:
if rand2 == -1 & rand3 == -1:
x2 = 5
y2 = 5
elif rand2 == -1 & rand3 > -1:
x2 = 10 - y1
elif rand2 > -1 & rand3 == -1:
y2 = 10 - x1
if x2<0:
x2 = 10 - y1
if y2<0:
y2 = 10 - x1
if x2>10:
x2 = 0
if y2>10:
y2 = 0
elif rand1 == 1:
x2 = x1
y2 = y1 + random.choice([-1, 1])
if x2 + y2 > 10:
y2 = 0
if y2 < 0:
y2 = 10 - x2
g1 = (4 ** x1 / math.factorial(x1)) * (4 ** y1 / math.factorial(y1))
g2 = (4 ** x2 / math.factorial(x2)) * (4 ** y2 / math.factorial(y2))
if g2 >= g1:
list_.append(np.array([x2, y2]))
elif g2 < g1:
u = np.random.rand()
if u < (g2 / g1):
list_.append(np.array([x2, y2]))
else:
list_.append(np.array([x1, y1]))
return list_

from this we got the following normal distribution:
![](https://i.imgur.com/ja8gPwC.png)

By comparing this to the distribution of the function we can see that this is an accurate representation of the expected distribution:
![](https://i.imgur.com/r2IusWL.png)

furthermore by using this method we improve our chisquared result.





## EXERCISE 7
_Nicholas Rose (s164580)_

Simulated annealing with the given functions for the travelling salesman problem

import numpy as np
import random
import seaborn as sns
import scipy.stats
import statistics
import math

def trav(m, k):
stations = []
s1 = 1
while len(stations) < len(m[:]):
stations.append(s1)
s1 += 1

stations.append(1)

cost = findcost(m, stations)

for i in range(k):
Tgay = T(k)

stations_proposal = permut(stations)

if findcost(m, stations_proposal) < findcost(m, stations):
stations = stations_proposal
else:
fprob = f(Tgay,m, stations, stations_proposal)
if random.uniform(0, 1) < fprob:
stations = stations_proposal
cost = findcost(m, stations)

return stations, cost

def T(k):
return 1/math.sqrt(1+k)

def f(T, m, S1, S2):
return math.exp(-(findcost(m,S2) - findcost(m,S1))/T)

def findcost(m, S): #m is matrix, S is proposed route
cost = 0
for i in range(len(S)-1):
cost = cost + m[(S[i]-1)][(S[i+1]-1)]
return cost

def permut(s): #switch two random elements of s

s_index = []
for i in range(len(s)-2):
s_index.append(i+1)

a = np.random.choice(s_index)
b = np.random.choice(s_index)

aa = s[a]
bb = s[b]

s[a] = bb
s[b] = aa

return s

## EXERCISE 8
_Nicholas Rose (s164580)_
To solve exercise 13 in Chapter 7 of Ross, we use the bootstrap approach to generate $r$ data sets by choosing at random from the given data set. Then, we find the mean for each of these data sets and then compute the variance of these means. Once we have a collection of means for each of the simulated data sets, we can see that these means fit a normal distribution when r is large:

![](https://i.imgur.com/50ngJXW.png)


We find the mean and standard deviation of these means and we can use this to find the probability that the underlying mean has a maximum difference of 5 to the observed mean. We use the built in normal distribution library to perform a two-tailed test. The reliability of our tests will increase with $r$. 

def boostrap(r):
data = [56, 101, 78, 67, 93, 87, 64, 72, 80, 69]
a = -5
b = 5
m = np.zeros((r, np.size(data)))
mu = np.zeros(r)
for i in range(r):
for j in range(np.size(data)):
m[i][j] = random.choice(data)

for i in range(r):
mu[i] = sum(m[i])/np.size(data)

mean = sum(mu)/r

var = statistics.variance(mu)

lower = scipy.stats.norm.cdf(mean-5, mean, math.sqrt(var))
prob = 1 - (2*lower)

return m, mu, mean, var

The following code defines the routine to return the bootstrap estimate of the variance of the median for a given data set. 

def boosts(data):

median_org = statistics.median(data)
mean_org = statistics.mean(data)

m = np.zeros((100, np.size(data)))
med = np.zeros(100)
meangu = np.zeros(100)
for i in range(100):
for j in range(np.size(data)):
m[i][j] = random.choice(data)

for i in range(100):
med[i] = statistics.median(m[i])
meangu[i] = statistics.mean(m[i])

meanmed = statistics.mean(med)
meanmean = statistics.mean(meangu)

varmed = statistics.variance(med)
varmean = statistics.variance(meangu)

return mean_org, median_org, meanmean, meanmed, varmean, varmed

Testing the method with a simulated Pareto distribution with N = 200 data values with underlying mean $\mu = 8.97940906077152$ and median $med = 0.9079021604107131$, after bootstrapping we find the following results:
$mean = 8.658400887314448$
$median = 0.9276452933918827$
$var(mean) =  20.328114005439403$
$var(median) = 0.013860234741911764$

As shown, the variance of the mean is significantly higher than the median, which suggests that the bootstrapping technique is more suited for estimating the median.