Department of Mathematics and Statistics MAST30028 Numerical Methods & Scientific Computing 2019
Assignment 1: Simulation & errors
Due: 5 p.m. september 5
This assignment is worth 20% of the total assessment in MAST30028.
Note: Submit copies of MATLAB programs and sufficient relevant output using the link on the LMS.
To aid marking, please start Questions 3 and 4 on a new page in your write-up.
1 Roulette [20 marks]
In Crown Casino a roulette wheel has 18 slots coloured red, 18 slots coloured black and 1 slot (0) coloured green. The red and black slots are also numbered from 1-36.
You can play various “games” or “systems” in roulette. Four possible games are:
• A. Betting on Red
This game involves just one bet. You bet $1 on red. If the ball lands on red you win $1 (and get your bet
back), otherwise you lose $1.
• B. Betting on a number
This game involves just one bet. You bet on a particular number, say 17; if the ball lands on that number
you win $35 ( and get your bet back), otherwise you lose $1.
• C. Martingale System
In this game you start by betting $1 on red. If you lose, you double your previous bet; if you win, you bet
$1 again. You continue to play until you have won $10, or the bet exceeds $100.
• D. Labouchere System
In this game you start with the list of numbers (1, 2, 3, 4). You bet the sum of the first and last numbers on red (initially $5). If you win you delete the first and last numbers from the list, which becomes (2,3), otherwise you add the sum to the end of your list (which becomes (1, 2, 3, 4, 5)). You repeat this process till your list is empty, or the bet exceeds $100. (If only one number is left on the list you bet that number).
Different games offer different playing experiences e.g. some allow you to win more often than you lose, some let you play longer, some cost more to play and some risk greater losses. Our aim is to compare the four game types above using:
a. the expected winnings/game
b. the proportion of games you win
c. the expected playing time (measured by the number of bets) d. the maximum amount you can lose
e. the maximum amount you can win
1
i.
ii.
iii. iv.
Write a MATLAB program which estimates (a), (b) and (c) for Games A-D by simulating 100,000 repe- titions. The program should play a different game, depending whether the input argument is ’A’ or ’B’ etc. The first line should be function roulette(gametype). There should be one subfunction for each type of game which simulates a single game of that type.
Check your program estimates for (a) and (b) for Game A by calculating the exact answers. What is the percentage error in your estimates for 100,000 repetitions, as measured by the mean of 5 runs, compared to the exact result?
Using mathematical arguments, calculate the exact answers for (d) and (e) for Games A-D. (Of course you might want to use your program to investigate this problem to help you find the theoretical answer).
Summarise your results in a table as follows:
In each entry of the first three columns of the table, enter the range of values (minimum, maximum) from 5 separate runs, each of 100,000 repetitions. Compare the variability of expected winnings in Games A and C, as measured by the half-range of 5 runs i.e. (max-min)/2. Compare the variability of expected winnings and expected playing time in Game D.
Buffon’s needle [10 marks]
Game
Exp. winnings
Prop wins
Exp. play time
Max loss
Max win
A
B
C
D
2
The following question was first considered by George Louis Leclerc, later Comte de Buffon, in 1733:
“If a thin, straight needle of length l is thrown at random onto the middle of a horizontal table ruled with parallel lines a distance d ≥ l apart, so that the needle lies entirely on the table, what is the probability that no line will be crossed by the needle?”
The answer depends on π−1 and so simulation of this experiment offers a way of estimating π−1. We will look at the complementary probability that the needle actually intersects with a ruled line on the table – this is called a crossing.
a. Theoretical answer We can think of the position of the needle as being determined by two random variables:
(i) Y = the perpendicular distance of the center of the needle from the nearest line on the table and
(ii) X = the angle that the top half of the needle makes with a ray through its centre, parallel to the
table lines and extending in a positive direction.
For the position of the needle to be random we require Y to be U(0, d) and X to be U(0,π). We then
2
define the sample space Ω of all possible outcomes or positions of the needle as Ω = [0, π] × [0, d ]. We
restrict attention to the case l ≤ d.
2
(i) Derive the inequality which X and Y must satisfy if the needle is to cross a ruled table line. Draw a picture of the sample space Ω and use your inequality to shade that part of it which corresponds to a crossing (will will refer to this region as the crossing region C).
(ii) As the needle is random the probability of falling in any region R in Ω can be calculated as the ratio of area of the region to the total area of Ω ie Area(R) . Using integration find the area of C and hence
dπ
2
confirm that the probability of a crossing is in fact 2l . πd
b. Simulation estimates Write a MATLAB function buffon.m to estimate the probability of a crossing directly using the proportion of N=100,000 simulated needle tosses which fall inside the crossing region C, for the case d = 2, l = 1. (You can use code similar to that used for Hit and miss Monte Carlo integration).
(i) Your function should compute the estimate the probability of a crossing, and a 2-sided 95% confidence interval for the probability. It should return a logical variable which is true if the exact value 1/π is contained in the CI and false otherwise.
(ii) By calling your function 100 times (from a driver function), check how many times the confidence interval does in fact contain the true probability. Is this as expected?
2
3 A bound on roundoff error [4 marks]
To analyse the roundoff error from n floating point multiplies, the following quantity arises: Πni=1(1 + δi) ≡ 1 + θn
where |δi| ≤ u and u is unit roundoff i.e. θn is the relative error after n floating point multiplies. Use mathematical induction to prove the bound
|θn| < nu 1−nu
assuming nu < 1. The same bound holds (with ≤) if we include divides as well, but you do not need to prove this.
4
a.
b. c.
d. e.
A second derivative.
Use Taylor series to derive the truncation error of the approximation f′′(x)≈ f(x+h)−2f(x)+f(x−h)
h2
assuming f ∈ C4.
Explain why dividing by h produces no roundoff error if h = 2−k,k ∈ Z.
[16 marks]
Assuming that built-in functions return exact answers but that addition/subtraction produce roundoff errors and that h = 2−k,k ∈ Z, show that the roundoff error RE in computing this expression satisfies the bound
|RE| ≤ (K1 + K2h)u/h2 where the constants K1,K2 depend on f and x.
Hence estimate, in terms of u, the optimal choice for h. You may assume hu ≪ u.
Write a script file called SecondDerivative.m to test your estimates for f(x) = exp(x) and x = 1,
h = 0.5, 0.25, . . . 2−52 and plot the errors, as in ForwardDifference.m. Comment on your results.
3