In [71]:
# Original NCL Code is often provided in comments.
# Some lines are directly modified to python syntax, but often it is written in more idiomatic python.
#
# ;***********************************************************
# ; resampling_1.ncl
# ;
# ; Concepts illustrated:
# ; - Generating random number seeds based on a 'clock'
# ; - Generate 'resampling with replacement' indices
# ; Use: generate_sample_indices( N, 1 )
# ; - Use histogram to display distribution
# ; - create a panel plot
# ;***********************************************************
# ;--- Specify sample sizes
# ;***********************************************************
NS = np.array([25, 100, 500, 1000, 5000, 10000, 100000, 1000000])
# NP = len(NS) # NCL --> not needed in Python version
# plot = new (NP, "graphic", "No_FillValue") #
#***********************************************************
#--- Generate and use random number seeds for initialization
#***********************************************************
# NCL defaults to fixed seeds for its random number generator,
# but that is not the case for Numpy, so re-seeding is not strictly necessary.
# Numpy has fairly advanced random number generation features (https://numpy.org/doc/stable/reference/random/index.html).
# NCL:
# rand1 = toint(systemfunc(" date +%s")) # --> NCL uses the system's date command
# Python:
# Method 1 (subprocess to use OS date just like NCL)
# import subprocess
# cmd = subprocess.run(['date', '+%s'], capture_output=True)
# rand1 = int(cmd.stdout)
# If you don't want to use the OS (directly) to get the date, you could use Python:
# Method 2 (datetime):
from datetime import datetime
curr_time = datetime.now()
rand1 = int(curr_time.strftime('%s'))
# Note: Method 1 and Method 2 produce the same result.
rand2 = int((54321 * rand1) % 2147483398)
# NCL : random_setallseed(rand1, rand2) # the ncl RNG uses two integers for seeding
# Python/Numpy:
# We can use a sequence (of any length) to seed.
# If no seed is given, "fresh, unpredictable entropy will be pulled from the OS" <-- Means we didn't need to worry about trying to make up a new seed.
rng = np.random.default_rng(seed=[rand1, rand2]) # rng = random number generator
# NCL Code from example
#***********************************************************
#--- Plot a histogram of index distributions for each sample size
#
#***********************************************************
# wks = gsn_open_wks ("png","resampling") ; send graphics to PNG file
# gsn_define_colormap(wks,"default") ; 30 distinct colors
# resh = True
# resh@gsnDraw = False
# resh@gsnFrame = False
# resh@tmXBLabelStride = 1
# resh@tiMainOffsetYF = -0.020 ; move tiMain down into plot
# ;resh@gsnHistogramBinWidth = ?? ; for 'tuning' ... if desired
#***********************************************************
#--- Loop
#***********************************************************
# do plotnum=0,2
# do np=0,NP-1
# iw := generate_sample_indices( NS(np), 1 )
# iw@long_name = "N="+NS(np)
# printMinMax(iw, 0)
# resh@tiMainString = iw@long_name
# delete(iw@long_name) ; do not want on plot
# if (np.eq.0 .or. np.eq.4) then ; reduce clutter
# resh@tiYAxisString = "Frequency"
# else
# resh@tiYAxisString = " " ; No title on Y axis.
# end if
# plot(np) = gsn_histogram(wks, iw ,resh)
# end do
# resP = True
# resP@gsnMaximize = True
# resP@gsnPanelMainString = "Distribution: Uniform Random Indices w/ Replacement"
# resP@gsnPanelLeft = 0.01 ; make room on left and right side of paneled
# resP@gsnPanelRight = 0.99 ; plots for text and tickmarks
# gsn_panel(wks,plot,(/2,4/),resP)
# end do
# Python using Numpy and Matplotlib
fig, ax = plt.subplots(figsize=(12,6), ncols=4, nrows=2, sharey=False, constrained_layout=True)
aa = ax.ravel() # makes a view of ax that can be indexed with one number.
iw = dict() # hold the results all together // This is optional, but could be useful if there are more steps than just making the histogram.
histograms = dict()
h_edges = dict()
for i, ns in enumerate(NS):
print(f"i = {i}, ns = {ns}")
iw[ns] = generate_sample_indices(ns, method=True, rng=rng)
histograms[ns], h_edges[ns] = np.histogram(iw[ns]) # keep histogram counts and bin edges.
aa[i].bar( h_edges[ns][0:-1] + np.diff(h_edges[ns]), histograms[ns], width=0.8*np.diff(h_edges[ns])[1] )
aa[i].set_title(f"N = {ns}")
plt.setp(ax[:,0], ylabel="FREQUENCY") # put labels on left side axes.
i = 0, ns = 25 i = 1, ns = 100 i = 2, ns = 500 i = 3, ns = 1000 i = 4, ns = 5000 i = 5, ns = 10000 i = 6, ns = 100000 i = 7, ns = 1000000
Out[71]:
[Text(0, 0.5, 'FREQUENCY'), Text(0, 0.5, 'FREQUENCY')]
