(* Student Name *)
(* July 24, 2018 *)
(* Physics III - Introduction to Statistical Physics *)
(* Assignment 5 - Problem 5)

(a) Mean and Standard Deviation of Binomial Distribution

For the binomial distribution where the probability of getting heads is 50%, the mean for N coin flips is < k> = N/2 and the standard deviation is σ_k= sqrt(N)/4. If we take N = 1000, we have <k> = 500 and σ_k = 15.811.

(b) Function Definitions

?For

?RandomReal

?If

?Print

?Range

?AppendTo

?Histogram

?Mean

?StandardDeviation

?N

(c) Simple Code Examples

Simple For-Loop

For[i = 0, i <5, i++, Print[i]](* Starting from i = 0, prints the number i and increments by 1. Repeats this process until i exceeds 5*)

Single Coin Toss Simulation

x=RandomReal[]; (* Selects a random real number between 0 and 1*)
If[x > 0.5,Print[H], Print[T]];(*If the randomly selected real number is greater than 0.5, then print H; otherwise print T *)

Three Coin Tosses Simulation

For[k=1,k<4,k++,x=RandomReal[]; (*Starting from k=1, selects a random real number, executes program and then increments by 1. Cycle repeates until k exceeds 4 *)
If[x > 0.5,Print[H], Print[T]]] (* If the randomly selected real number is greater than 0.5, then print H; otherwise print T*)

Element Extraction

vec = {2,4,6,8,10,12,14,16,18,20};  (* Defines a vector of even numbers from 0 to 20*)
vec[[4]](* Selects the fourth element from the vector of numbers *)

Element appending

vec = {1, 2, 3}; (* Defines a vector consisting of the numbers 1, 2, 3 in sequence *)
Print[vec] (* Prints the vector *)
AppendTo[vec, 4]; (* Appends to vector the number 4 and thus extends the length by 1 *)
Print[vec] (* Print the new vector *)

(Simulating Coin Flips)

(d) Single Coin-Flip Simulations

Nflips = 1000; (*Student Mistake in Assignment;  Defines the total number of coin flips for a single simulation *)
flip = Range[0,Nflips-1]; (*Defines numbered sequence of flips from 0 to Nflips-1 *)
HeadsCount ={0}; (*Sets the initial count of heads, for zero coin flips to be zero*)
For[k=1,k<Nflips,k++,
x=RandomReal[]; If[x > 0.5,{AppendTo[HeadsCount ,HeadsCount [[k]]+1]}, 
{AppendTo[HeadsCount ,HeadsCount [[k]] ]}]] (*Student Mistake in Assignment;  gets a random real number; if number is greater than 0.5, make the next count of number of heads, "1+previous count"; otherwise make next count of number of heads the same as previous count; repeats this procedure a total of Nflips times, and thereby obtains a count of heads obtained as a function of current flip number  *)
ListLinePlot[Transpose[{flip,HeadsCount }]] (*Plots the count of total heads as a function of the current flip*)

(e) Many Coin-Flip Simulations

Nsimulations=5000; (*Student Mistake in Assignment; Defines the number of simulations, that is the number of times we perform 1000 coin flips*)
Nflips=1000; (*Student Mistake in Assignment;  Defines the number of coin flips in a single simulation. *)
flip = Range[0,Nflips-1]; (* Defines numbered sequence of flips from 0 to Nflips-1 *)
SimHeadsCount={}; (* Empty list which will contain the total count of heads for each simulation *)
For[j=1,j<Nsimulations,j++,HeadsCount ={0};For[k=1,k<Nflips,k++,x=RandomReal[];If[x > 0.5,AppendTo[HeadsCount,HeadsCount[[k]]+1],AppendTo[HeadsCount,HeadsCount[[k]] ]]];AppendTo[SimHeadsCount,HeadsCount[[Nflips] ]]]; (* Student Mistake in Assignment;  gets a random real number; if number is greater than 0.5, make the next count of number of heads, "1+previous count"; otherwise make next count of number of heads the same as previous count; repeats this procedure a total of Nflips times, and thereby obtains a count of heads obtained as a function of current flip number. We then repeat this procedure Nflips times *)

Histogram[SimHeadsCount,50] (* Creates a histogram of the total number of heads obtained across all simulations *)

N[Mean[SimHeadsCount]] (* Computes the mean of the distribution *)

N[StandardDeviation[SimHeadsCount]](* Computes the standard deviation of the distribution *)

(f) Describing the Distribution

(The distribution looks like a Gaussian distribution. By the computations at the end of the simulation, we find, 
mean = 499.301
standard deviation = 15.922

In part (a) we found (or quoted) the theoretical result <k>=500 and σ_k=15.811 for a Binomial distribution with p=1/2 and N=1000. Therefore the results of the simulation are close to the theoretical results. We can therefore consider the distribution of heads for 1000 coin flips to be a Gaussian distribution with mean 500 and standard deviation 15.811)