N = np .
array ( N )
Statistical Mechanics Practical Solutions def model (t , N0 , lam ) :
return N0 * np . exp ( - lam * t )
popt , _ = curve_fit ( model , t , N )
plt . scatter (t , N )
plt . plot (t , model (t , * popt ) )
Question 1(a): Uniform Random Sample plt . show ()
Generate 200 uniform random numbers, scatter plot, correlogram. print ( popt )
import numpy as np
import matplotlib . pyplot as plt
Question 2(a): Uniform Histogram
x = np . random . uniform (0 , 1 , 200)
import numpy as np
plt . scatter ( x [:100] , x [100:])
import matplotlib . pyplot as plt
plt . xlabel ( " Subsample 1 " )
plt . ylabel ( " Subsample 2 " )
X = np . random . uniform (0 , 1 , 10000)
plt . show ()
plt . hist (X , bins =20 , density = True )
lags = 20
plt . show ()
acf = [ np . corrcoef ( x [: - k ] , x [ k :]) [0 ,1] for k in range (1 , lags ) ]
Y = - np . log (1 - X )
plt . stem ( range (1 , lags ) , acf )
plt . show ()
plt . hist (Y , bins =20 , density = True )
plt . show ()
Question 1(b): Radioactive Decay
Decay law:
Question 2(b): 2D Random Walk
ln 2 RRMS ∝ t1/2
N (t) = N0 e−λt , λ=
500
import numpy as np
import numpy as np import matplotlib . pyplot as plt
import matplotlib . pyplot as plt from scipy . optimize import curve_fit
from scipy . optimize import curve_fit
steps = 1000
N0 = 10000 trials = 500
half_life = 500 R = []
lam = np . log (2) / half_life
for t in range (1 , steps ) :
t = np . arange (0 , 3000 , 10) r2 = []
N = [ N0 ] for _ in range ( trials ) :
x = np . random . choice ([ -1 ,1] , t )
for _ in t [1:]: y = np . random . choice ([ -1 ,1] , t )
decay = np . random . binomial ( N [ -1] , lam *10) r2 . append (( x . sum () ) **2 + ( y . sum () ) **2)
N . append ( N [ -1] - decay ) R . append ( np . sqrt ( np . mean ( r2 ) ) )
1 2
Question 4(a): Fermi-Dirac Distribution
R = np . array ( R )
time = np . arange (1 , steps ) 1
f (E) =
e(E−EF )/kT + 1
def power (t , A , n ) :
import numpy as np
return A * t ** n
import matplotlib . pyplot as plt
popt , _ = curve_fit ( power , time [: len ( R ) ] , R )
k = 8.617 e -5
Ef = 0.15
plt . scatter ( time [: len ( R ) ] , R )
plt . plot ( time [: len ( R ) ] , power ( time [: len ( R ) ] , * popt ) )
x = np . arange ( -3 , 3 , 0.1)
plt . show ()
E = Ef + x
print ( popt )
def f (E , T ) :
return 1 / ( np . exp (( E - Ef ) /( k * T ) ) + 1)
Question 3(a): Coin Toss (6 Tosses) plt . plot (x , f (E , 3) , label = " 3 K " )
plt . plot (x , f (E , 300) , label = " 300 K " )
6 plt . legend ()
2
P (2) = plt . show ()
64
import numpy as np
import matplotlib . pyplot as plt
from math import comb Question 4(b): Monte Carlo Integration
heads = np . random . binomial (6 , 0.5 , 10000) import numpy as np
plt . hist ( heads , bins =7 , density = True ) N = 100000
plt . show () a, b = 0, 2
P_exp = np . mean ( heads == 2) x = np . random . uniform (a , b , N )
P_theory = comb (6 ,2) / 64 f = x **2
print ( P_exp , P_theory ) I = (b - a ) * np . mean ( f )
print ( I )
Question 3(b): Exponential Moments
Question 5(a): Standard Normal
import numpy as np
import numpy as np
X = np . random . uniform (0 ,1 ,1000)
Y = - np . log (1 - X ) x = np . random . normal (0 ,1 ,10000)
mean = np . mean ( Y ) print ( np . mean ( x ) )
mu2 = np . mean (( Y - mean ) **2) print ( np . var ( x ) )
mu3 = np . mean (( Y - mean ) **3)
print ( np . mean ( np . abs ( x ) <=1) )
print ( mu2 , mu3 ) print ( np . mean ( np . abs ( x ) <=3) )
3 4
Question 5(b): Coin Flip Game
import numpy as np
wins = 0
M = 10000
for _ in range ( M ) :
while True :
if np . random . rand () < 0.5:
wins += 1
break
if np . random . rand () < 0.5:
break
print ( wins / M )
Question 6(b): CLT Normal Distribution
import numpy as np
import matplotlib . pyplot as plt
X = np . random . uniform (0 ,1 ,(10000 ,12) )
Z = X . sum ( axis =1) - 6
plt . hist (Z , bins =30 , density = True )
plt . show ()
print ( np . mean ( Z ) , np . var ( Z ) )