Working with Data from Textbook¶
Data sets can be treated in two ways in python:
- Numpy arrays
- Pandas DataFrame
Numpy arrays are vectors and matrices. If data from the textbook are small in size, importing them to Numpy arrays works well. Statistics can computed using
Pandas DataFrames are similar to tabulation of data in spreadsheet files, with the first row or the first several rows as the data header (variable names), and a column represents the values of a variable. Each row is one observation.
Whenever it is sufficient, we use Numpy and Scipy.stats for computing statistics.
Note: Numpy, Scipy, and Pandas are not standard libraries of Python, hence must be installed separately and must be loaded into the Python activation.
Reading data as Pandas DataFrame¶
In [41]:
import numpy as np
import pandas as pd
pd.set_option("mode.copy_on_write", True)
filepath = "/home/xcst/teaching/Practice_Statistics/dataSets/psbe5e_csv/Chapter1/"
tbill=pd.read_csv(filepath+"eg01-12tbill.csv", index_col=0,date_format="%m/%d/%Y")
tbill.index
Out[41]:
DatetimeIndex(['1958-12-12', '1958-12-19', '1958-12-26', '1959-01-02',
'1959-01-09', '1959-01-16', '1959-01-23', '1959-01-30',
'1959-02-06', '1959-02-13',
...
'2018-01-26', '2018-02-02', '2018-02-09', '2018-02-16',
'2018-02-23', '2018-03-02', '2018-03-09', '2018-03-16',
'2018-03-23', '2018-03-30'],
dtype='datetime64[ns]', name='Date', length=3095, freq=None)
In [42]:
tbill.describe()
Out[42]:
| Rate | |
|---|---|
| count | 3095.000000 |
| mean | 4.730753 |
| std | 3.135163 |
| min | 0.030000 |
| 25% | 2.745000 |
| 50% | 4.830000 |
| 75% | 6.360000 |
| max | 15.760000 |
Plot histogram of t-bill rates¶
In [43]:
import matplotlib.pyplot as plt
import matplotlib as mpl
# Customize plot style sheet, rc stands for runtime configuration
mpl.rc('axes', labelsize=16)
mpl.rc('xtick', labelsize=16)
mpl.rc('ytick', labelsize=16)
plt.rc('font', size=18)
plt.rc('axes', labelsize=16, titlesize=20)
plt.rc('legend', fontsize=18)
plt.rc('xtick', labelsize=16)
plt.rc('ytick', labelsize=16)
plt.rcParams['lines.linewidth'] = 2.5
#plt.rcParams['axes.titlesize'] = 20
plt.rcParams["figure.figsize"] = (16/1.5, 9/1.5) #set default figure size
In [44]:
fig, ax = plt.subplots(nrows=1, ncols=1, sharey='all')
ax.hist(tbill['Rate'], bins=16, linewidth=0.5, edgecolor="white")
plt.autoscale(enable=True, axis='x', tight=True)
ax.set_title('Histogram of T-bill interest rates')
ax.tick_params(axis='x')
ax.tick_params(axis='y')
ax.set_xlabel("Rate")
ax.set_ylabel('Count')
#ax.legend(loc='upper right', fontsize=18)
fig.tight_layout()
filename = 'tbill_histogram.png'
plt.savefig(filename, dpi=300, format='png')
plt.show()
Plot the time series of T-bill interest rates¶
In [45]:
# time series of T-bill interest rate
starttime = '1970-01-02'
fig, ax = plt.subplots(nrows=1, ncols=1, sharey='all')
ax.plot(tbill[starttime:]['Rate'],color="green", linestyle="-",label='T-bill Rate')
plt.autoscale(enable=True, axis='x', tight=True)
ax.set_title('Historical T-bill Interest Rates')
ax.tick_params(axis='x')
ax.tick_params(axis='y')
ax.set_xlabel("Year")
ax.set_ylabel('T-bill Rate')
ax.legend(loc='upper right', fontsize=18)
fig.tight_layout()
filename = 'tbill_history.png'
plt.savefig(filename, dpi=300, format='png')
plt.show()
More example: the Labor Force Survey¶
We use the public use micro files (PUMFs) of the Labour Force Survey (LFS) 2024 May data. LFS is a monthly survey on employment status in Canada conducted by Statistics Canada, it is the main source of statistics on the Canadian labor market condition. For example, the unemployment rate in headline news is calculated from the LFS.
In the LFS, an observation (individual) has a weight, indicating how many Canadians are closely similar to the individual surveyed by the LFS. For the purpose of this course, we ignore the weight, which means any statistics calculated here are likely different from those reported by Statistics Canada.
For variables included the LFS data file, please refer to file LFS_PUMF_EPA_FGMD_codebook.csv.
Import libraries¶
We already imported libraries Numpy, Pandas, and Matplotlib in the above. We do it again here.
In [46]:
import numpy as np
import pandas as pd
pd.set_option("mode.copy_on_write", True)
In [47]:
import matplotlib.pyplot as plt
import matplotlib as mpl
# Customize plot style sheet, rc stands for runtime configuration
mpl.rc('axes', labelsize=16)
mpl.rc('xtick', labelsize=16)
mpl.rc('ytick', labelsize=16)
plt.rc('font', size=18)
plt.rc('axes', labelsize=16, titlesize=20)
plt.rc('legend', fontsize=18)
plt.rc('xtick', labelsize=16)
plt.rc('ytick', labelsize=16)
plt.rcParams['lines.linewidth'] = 2.5
#plt.rcParams['axes.titlesize'] = 20
plt.rcParams["figure.figsize"] = (16/1.5, 9/1.5) #set default figure size
Read in LFS 2024 May¶
In [48]:
# data file is in Stata format, created by Shutao
lfs=pd.read_stata("./datafiles/LFS202405.dta")
lfs.columns
Out[48]:
Index(['rec_num', 'year', 'month', 'lfsstat', 'province', 'cma', 'agegrp',
'sex', 'marstat', 'educ', 'mjh', 'ftptlast', 'cowmain', 'immig',
'naics43', 'soc80_49', 'soc80_21', 'uhrsmain', 'ahrsmain', 'ftptmain',
'uhrsall', 'ahrsall', 'xtrahrs', 'whyptold', 'whyptnew', 'tenure',
'prevten', 'hrlyearn', 'union', 'permtemp', 'estsize', 'firmsize',
'durunemp', 'flowunem', 'durjless', 'availabl', 'schooln', 'relrefn',
'efamtype', 'efamsize', 'efamempl', 'efamunem', 'sp_age', 'sp_lfsst',
'sp_uhrsm', 'sp_uhrst', 'sp_cowm', 'agyownkn', 'finalwt', 'naics21',
'naics18', 'nocs2006_25', 'nocs2006_47', 'nocs2016_10', 'nocs2016_40',
'empstat', 'empclass', 'laborforce', 'employed', 'unemployed',
'busSector'],
dtype='object')
Describing categorical variables¶
Use age group as an example. From codebook, we note that there are 12 age groups.
In [49]:
lfs.agegrp.unique()
Out[49]:
array([12., 4., 1., 11., 5., 7., 6., 8., 2., 3., 10., 9.],
dtype=float32)
In [50]:
count = lfs['agegrp'].value_counts(sort=True, ascending=True, dropna=True)
count2 = count.sort_index()
print(count2)
shares = 100*count2/count2.sum()
shares
agegrp 1.0 7561 2.0 7327 3.0 7637 4.0 8730 5.0 9056 6.0 9038 7.0 8081 8.0 8293 9.0 8515 10.0 9953 11.0 9046 12.0 19734 Name: count, dtype: int64
Out[50]:
agegrp 1.0 6.692868 2.0 6.485735 3.0 6.760142 4.0 7.727647 5.0 8.016217 6.0 8.000283 7.0 7.153163 8.0 7.340822 9.0 7.537333 10.0 8.810226 11.0 8.007365 12.0 17.468200 Name: count, dtype: float64
In [51]:
xlabels=['15-19', '20-24', '25-29','30-34','35-39','40-44',
'45-49','50-54','55-59','60-64','65-69','70+']
x = np.arange(len(xlabels)) # the label locations
width = 0.75 # the width of the bars
fig, ax = plt.subplots(nrows=1, ncols=1, sharey='all')
rects1 = ax.bar(x, shares.round(1), width)
# Add some text for labels, title and custom x-axis tick labels, etc.
ax.set_title('Age distribution working-age population, 2024 May')
ax.set_ylabel('Percent')
plt.xlabel("Data source: Labour Force Survey PUMFs")
ax.set_xticks(x, xlabels)
ax.tick_params(axis='x')
#ax.legend()
# ax.bar_label(rects1, padding=2)
fig.tight_layout()
filename = 'age_barShare.png'
plt.savefig(filename, dpi=400, format='png')
plt.show()
In [52]:
xlabels=['15-19', '20-24', '25-29','30-34','35-39','40-44',
'45-49','50-54','55-59','60-64','65-69','70+']
x = np.arange(len(xlabels)) # the label locations
width = 0.75 # the width of the bars
fig, ax = plt.subplots(nrows=1, ncols=1, sharey='all')
ax.set_title('Age distribution working-age population, 2024 May')
ax.pie(shares, labels=xlabels, autopct='%1.1f%%', radius=1.0)
fig.tight_layout()
filename = 'age_pieShare.png'
plt.savefig(filename, dpi=400, format='png')
plt.show()
Describing quantitative variables¶
We use the variable hrlyearn as an example, it is the hourly earnings for individuals that are employees. Self-employed individuals are excluded as their income is not salary.
In [53]:
lfs.hrlyearn.describe()
Out[53]:
count 58377.000000 mean 34.373962 std 18.625530 min 6.870000 25% 21.000000 50% 29.230000 75% 43.000000 max 205.130005 Name: hrlyearn, dtype: float64
In [54]:
fig, ax = plt.subplots(nrows=1, ncols=1, sharey='all')
ax.hist(lfs['hrlyearn'], bins=15, linewidth=0.5, edgecolor="white")
plt.autoscale(enable=True, axis='x', tight=True)
ax.set_title('Distribution of Hourly Earnings, 2024 May')
ax.tick_params(axis='x')
ax.tick_params(axis='y')
plt.xlabel("\nData source: Labour Force Survey PUMFs")
ax.set_ylabel('Number of individuals')
#ax.legend(loc='upper right', fontsize=18)
fig.tight_layout()
filename = 'hourlyEarnings.png'
plt.savefig(filename, dpi=300, format='png')
plt.show()
In [55]:
fig, ax = plt.subplots(nrows=1, ncols=1, sharey='all')
ax.hist(lfs['hrlyearn'], bins=15, density=True,linewidth=0.5, edgecolor="white")
plt.autoscale(enable=True, axis='x', tight=True)
ax.set_title('Distribution of Hourly Earnings, 2024 May')
ax.tick_params(axis='x')
ax.tick_params(axis='y')
plt.xlabel("\nData source: Labour Force Survey PUMFs")
ax.set_ylabel('Fraction of individuals')
#ax.legend(loc='upper right', fontsize=18)
fig.tight_layout()
filename = 'hourlyEarnings.png'
plt.savefig(filename, dpi=300, format='png')
plt.show()
Time series¶
We plot the unemployment rate by sex in Canada. The data source is Statistics Canada Table 14100287, which is a major statistics table produced from the Labour Force Survey (LFS).
I already downloaded the raw data file and processed it.
In [56]:
#---read in LFS table: U, E, P
import pandas as pd
filename = './datafiles/tbl14100287Final3.csv'
lfs2 = pd.read_csv(filename, index_col=0)
lfs2.index = pd.to_datetime(lfs2.index)
lfsQtr = lfs2.resample('QS-OCT').mean()
lfsQtr.columns
Out[56]:
Index(['emplBothSex', 'emplFemale', 'emplFullBothSex', 'emplFullFemale',
'emplFullMale', 'emplMale', 'emplPartBothSex', 'emplPartFemale',
'emplPartMale', 'emplRateBothSex', 'emplRateFemale', 'emplRateMale',
'labForceBothSex', 'labForceFemale', 'labForceMale',
'participRateBothSex', 'participRateFemale', 'participRateMale',
'populationBothSex', 'populationFemale', 'populationMale',
'unemplBothSex', 'unemplFemale', 'unemplMale', 'unemplRateBothSex',
'unemplRateFemale', 'unemplRateMale'],
dtype='object')
In [57]:
lfsQtr.head()
Out[57]:
| emplBothSex | emplFemale | emplFullBothSex | emplFullFemale | emplFullMale | emplMale | emplPartBothSex | emplPartFemale | emplPartMale | emplRateBothSex | ... | participRateMale | populationBothSex | populationFemale | populationMale | unemplBothSex | unemplFemale | unemplMale | unemplRateBothSex | unemplRateFemale | unemplRateMale | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| ref_date | |||||||||||||||||||||
| 1976-01-01 | 9666.900000 | 3566.333333 | 8469.300000 | 2735.100000 | 5734.200000 | 6100.533333 | 1197.566667 | 831.233333 | 366.300000 | 57.233333 | ... | 77.933333 | 16891.700000 | 8538.533333 | 8353.200000 | 718.166667 | 307.766667 | 410.400000 | 6.933333 | 7.933333 | 6.300000 |
| 1976-04-01 | 9737.533333 | 3604.933333 | 8529.633333 | 2756.800000 | 5772.833333 | 6132.600000 | 1207.866667 | 848.100000 | 359.733333 | 57.266667 | ... | 77.733333 | 17008.100000 | 8599.700000 | 8408.400000 | 718.133333 | 312.233333 | 405.900000 | 6.866667 | 7.966667 | 6.200000 |
| 1976-07-01 | 9778.400000 | 3635.066667 | 8562.966667 | 2779.300000 | 5783.700000 | 6143.366667 | 1215.433333 | 855.833333 | 359.666667 | 57.100000 | ... | 77.533333 | 17121.433333 | 8659.100000 | 8462.333333 | 753.400000 | 336.066667 | 417.333333 | 7.166667 | 8.433333 | 6.366667 |
| 1976-10-01 | 9824.700000 | 3673.900000 | 8568.266667 | 2790.333333 | 5777.900000 | 6150.800000 | 1256.433333 | 883.600000 | 372.866667 | 57.066667 | ... | 77.533333 | 17210.966667 | 8705.666667 | 8505.333333 | 785.966667 | 338.466667 | 447.500000 | 7.433333 | 8.433333 | 6.766667 |
| 1977-01-01 | 9869.733333 | 3700.566667 | 8605.300000 | 2810.800000 | 5794.500000 | 6169.166667 | 1264.433333 | 889.733333 | 374.666667 | 57.033333 | ... | 77.700000 | 17300.700000 | 8752.766667 | 8548.000000 | 829.933333 | 358.833333 | 471.100000 | 7.766667 | 8.833333 | 7.100000 |
5 rows × 27 columns
In [58]:
#------unemployment rate by sex
startt = '2000-01-01'
fig, ax = plt.subplots(nrows=1, ncols=1, sharey='all')
ax.plot(lfsQtr[startt:]['unemplRateFemale'],color="blue", linestyle="-",label='Female')
ax.plot(lfsQtr[startt:]['unemplRateMale'],color="red", linestyle="--",label='Male')
plt.autoscale(enable=True, axis='x', tight=True)
ax.set_title('Unemployment Rate by Sex in Canada')
ax.tick_params(axis='x')
ax.tick_params(axis='y')
plt.xlabel("\nData source: Statistics Canada Table 14100287")
ax.set_ylabel('Percentage')
ax.legend(loc='upper left', fontsize=18)
fig.tight_layout()
filename = 'unempBySex_Canada.png'
plt.savefig(filename, dpi=300, format='png')
plt.show()
Density curves¶
We plot the density curve of hourly earnings in the LFS.
In [59]:
import seaborn as sns
fig, ax = plt.subplots(nrows=1, ncols=1, sharey='all')
sns.histplot(lfs['hrlyearn'], bins=40, stat='density',kde=True, color='blue', edgecolor='red')
plt.autoscale(enable=True, axis='x', tight=True)
ax.set_title('Distribution of Hourly Earnings, 2024 May')
ax.tick_params(axis='x')
ax.tick_params(axis='y')
plt.xlabel("\nData source: Labour Force Survey PUMFs")
ax.set_ylabel('Fraction of individuals')
#ax.legend(loc='upper right', fontsize=18)
fig.tight_layout()
filename = 'hourlyEarnings_density.png'
plt.savefig(filename, dpi=300, format='png')
plt.show()
In [60]:
lfs['log_hrlyearn'] = np.log(lfs['hrlyearn'])
fig, ax = plt.subplots(nrows=1, ncols=1, sharey='all')
sns.histplot(lfs['log_hrlyearn'], bins=20, stat='density',kde=True, color='blue', edgecolor='red')
plt.autoscale(enable=True, axis='x', tight=True)
ax.set_title('Distribution of Logaritham of Hourly Earnings, 2024 May')
ax.tick_params(axis='x')
ax.tick_params(axis='y')
plt.xlabel("\nData source: Labour Force Survey PUMFs")
ax.set_ylabel('Fraction of individuals')
#ax.legend(loc='upper right', fontsize=18)
fig.tight_layout()
filename = 'hourlyEarnings_logDensity.png'
plt.savefig(filename, dpi=300, format='png')
plt.show()
Normal distribution¶
We calculate cumulative proportion of SAT scores that are at least 820, following example on slides.
In [61]:
from scipy.stats import norm
mu = 1029.0
sigma = 209.0
print("Proportion of students whose SAT scores are at least 820 is: ", 1.0 - norm.cdf(820,mu,sigma))
Proportion of students whose SAT scores are at least 820 is: 0.8413447460685429
In [62]:
z = (820.0-mu)/sigma
print(f'The z-score of an SAT score of 820 is {z:.2f}')
print(f"Proportion of students for whom the z-score of SAT scores is at least {z} is:", 1.0 - norm.cdf(z))
The z-score of an SAT score of 820 is -1.00 Proportion of students for whom the z-score of SAT scores is at least -1.0 is: 0.8413447460685429
Normal Quantile Plot¶
Normal quantile plot is also called Q-Q plot. It's a visual tool to decide how close a distribution is to normal distribution.
Idea: If a variable has a normal distribution, its z-score at each percentile should be equal the percentile of standard normal distribution.
Import IQ data as Numpy array¶
In [63]:
# Import IQ data as Numpy array
filepath = "/home/xcst/teaching/Practice_Statistics/dataSets/psbe5e_csv/Chapter1/"
IQ = np.genfromtxt(filepath+'eg01-19iq.csv',skip_header=1)
Step 1. Sort values from smallest to the largest value¶
In [64]:
IQ_sorted = np.sort(IQ)
print("sorted IQ",np.round(IQ_sorted,1))
sorted IQ [ 81. 82. 89. 90. 94. 96. 97. 100. 101. 101. 101. 102. 102. 102. 103. 105. 106. 108. 109. 109. 109. 110. 110. 110. 112. 112. 113. 113. 114. 114. 114. 115. 116. 117. 117. 117. 118. 118. 122. 122. 123. 124. 124. 124. 125. 126. 127. 127. 128. 130. 131. 133. 134. 134. 136. 137. 139. 139. 142. 145.]
In [65]:
np.size(IQ_sorted)
Out[65]:
60
In [66]:
np.shape(IQ_sorted)
Out[66]:
(60,)
Step 2. Obtain the percentage of each sorted IQ value¶
Each sorted IQ value is a percentile, corresponding to it is the percent. For example, for IQ=96, it is the 6th value. 6/60=0.1 is the percentage of percentile value 96, meaning that 10% of IQ values are at or below 96. That is, 96 is the 10th percentile.
In [67]:
from scipy import stats
IQ_pct = stats.mstats.plotting_positions(IQ_sorted, alpha=0, beta=1)
# IQ_pct = stats.mstats.plotting_positions(IQ_sorted, alpha=0.4, beta=0.4)
print("Proportion of IQ values equal to or below each IQ value:")
print(np.round(IQ_pct,3))
Proportion of IQ values equal to or below each IQ value: [0.017 0.033 0.05 0.067 0.083 0.1 0.117 0.133 0.15 0.167 0.183 0.2 0.217 0.233 0.25 0.267 0.283 0.3 0.317 0.333 0.35 0.367 0.383 0.4 0.417 0.433 0.45 0.467 0.483 0.5 0.517 0.533 0.55 0.567 0.583 0.6 0.617 0.633 0.65 0.667 0.683 0.7 0.717 0.733 0.75 0.767 0.783 0.8 0.817 0.833 0.85 0.867 0.883 0.9 0.917 0.933 0.95 0.967 0.983 1. ]
Step 3. For each proportion (the percentage) displayed above, find the percentile according to the standard normal distribution.¶
These percentiles are z-scores obtained from standard normal distribution.
For example, we want to find the 10th percentile on the standard normal distribution.
In [68]:
normal_scores = stats.norm.ppf(IQ_pct)
print("Normal z-scores corresponding to the proportions of IQ values:")
print(np.round(normal_scores,4))
Normal z-scores corresponding to the proportions of IQ values: [-2.128 -1.8339 -1.6449 -1.5011 -1.383 -1.2816 -1.1918 -1.1108 -1.0364 -0.9674 -0.9027 -0.8416 -0.7835 -0.7279 -0.6745 -0.6229 -0.573 -0.5244 -0.477 -0.4307 -0.3853 -0.3407 -0.2967 -0.2533 -0.2104 -0.1679 -0.1257 -0.0837 -0.0418 0. 0.0418 0.0837 0.1257 0.1679 0.2104 0.2533 0.2967 0.3407 0.3853 0.4307 0.477 0.5244 0.573 0.6229 0.6745 0.7279 0.7835 0.8416 0.9027 0.9674 1.0364 1.1108 1.1918 1.2816 1.383 1.5011 1.6449 1.8339 2.128 inf]
Step 4. Plot normal scores (theoretical) against the sorted IQ data¶
If IQ follows a normal distribution, we expect a straight line.
In [69]:
area = 0.10*(IQ.max())
fig, ax = plt.subplots(nrows=1, ncols=1,sharey='all')
ax.scatter(normal_scores, IQ_sorted, s=area, marker='+', c='r',alpha=0.8, label='Data')
# ax.plot(normal_scores, normal_scores, color="red", linewidth=2.5, linestyle="-",label='Theoretical')
plt.autoscale(enable=True, axis='x', tight=False)
ax.set_title('Normal Quantile Plot for IQ', fontsize=20)
ax.tick_params(axis='x', labelsize=16)
ax.tick_params(axis='y', labelsize=16)
plt.xlabel("Normal Scores", fontsize=16)
ax.set_ylabel('IQ', fontsize=16)
ax.legend(loc='upper left', fontsize=16)
#fig.tight_layout()
plt.show()
plt.savefig('NQplot_IQ.png',dpi=300, format='png')
<Figure size 1066.67x600 with 0 Axes>
It's approximately a straight line, so IQ has approximately a normal distribution.
Note: We could have used the z-scores of IQ values on the vertical axis, then the plotted line will be a 45-degree line if IQs follow a normal distribution.
Step 4 alternative. Plot normal scores (theoretical) against the z-scores of sorted IQ data¶
If IQ follows a normal distribution, we expect this plot is a 45-degree line.
To obtain the z-score of IQ data, we need estimates of standard deviation and the mean value.
In [70]:
# Mean and standard deviation of IQ data
IQ_mean = np.mean(IQ_sorted)
IQ_std = np.sqrt(sum((IQ_sorted - IQ_mean)**2)/(np.size(IQ_sorted)-1))
IQ_zscore = (IQ_sorted - IQ_mean) / IQ_std
IQ_zscore
Out[70]:
array([-2.29602697e+00, -2.22846365e+00, -1.75552038e+00, -1.68795705e+00,
-1.41770376e+00, -1.28257711e+00, -1.21501378e+00, -1.01232381e+00,
-9.44760485e-01, -9.44760485e-01, -9.44760485e-01, -8.77197161e-01,
-8.77197161e-01, -8.77197161e-01, -8.09633836e-01, -6.74507188e-01,
-6.06943863e-01, -4.71817215e-01, -4.04253890e-01, -4.04253890e-01,
-4.04253890e-01, -3.36690566e-01, -3.36690566e-01, -3.36690566e-01,
-2.01563918e-01, -2.01563918e-01, -1.34000593e-01, -1.34000593e-01,
-6.64372689e-02, -6.64372689e-02, -6.64372689e-02, 1.12605541e-03,
6.86893797e-02, 1.36252704e-01, 1.36252704e-01, 1.36252704e-01,
2.03816028e-01, 2.03816028e-01, 4.74069326e-01, 4.74069326e-01,
5.41632650e-01, 6.09195974e-01, 6.09195974e-01, 6.09195974e-01,
6.76759298e-01, 7.44322623e-01, 8.11885947e-01, 8.11885947e-01,
8.79449271e-01, 1.01457592e+00, 1.08213924e+00, 1.21726589e+00,
1.28482922e+00, 1.28482922e+00, 1.41995587e+00, 1.48751919e+00,
1.62264584e+00, 1.62264584e+00, 1.82533581e+00, 2.02802578e+00])
In [71]:
plt.rcParams["figure.figsize"] = (10,8) #re-set the figure dimension size
area = 15*(IQ_zscore.max())
fig, ax = plt.subplots(nrows=1, ncols=1,sharey='all')
ax.scatter(normal_scores, IQ_zscore, s=area, marker='+', c='r',alpha=0.8,label="IQ-zscore")
ax.plot(normal_scores, normal_scores, color="navy", linestyle="-",label='45-degree line')
plt.autoscale(enable=True, axis='x', tight=False)
ax.set_title('Normal Quantile Plot for IQ', fontsize=20)
ax.tick_params(axis='x', labelsize=16)
ax.tick_params(axis='y', labelsize=16)
plt.xlabel("Normal z-scores", fontsize=16)
ax.set_ylabel('IQ z-scores', fontsize=16)
# ax.legend(loc='upper left', fontsize=16)
#fig.tight_layout()
plt.show()
plt.savefig('NQplot_IQzscore.png',dpi=300, format='png')
<Figure size 1000x800 with 0 Axes>
More example on Normal Quantile Plot¶
In [72]:
import numpy as np
import pandas as pd
pd.set_option("mode.copy_on_write", True)
filepath = "/home/xcst/teaching/Practice_Statistics/dataSets/psbe5e_csv/Chapter1/"
df=pd.read_csv(filepath+"ex01-56stout.csv", index_col=0,date_format="%m/%d/%Y")
df.index
Out[72]:
Index(['Bristol', 'Cardiff', 'English Agents', 'English O', 'English P',
'English R', 'Glasgow', 'Liverpool', 'London', 'Manchester',
'Newcastle-on-Tyne', 'Scottish'],
dtype='object', name='Region')
In [73]:
import numpy as np
a = 3.
rng = np.random.default_rng()
s = rng.pareto(a, 10)
s
Out[73]:
array([0.00458193, 0.14230026, 0.33959405, 0.24348738, 0.39992794,
0.02721935, 0.02276702, 0.07797975, 0.96993578, 0.77630536])
In [74]:
s_sorted = np.sort(s)
s_sorted
Out[74]:
array([0.00458193, 0.02276702, 0.02721935, 0.07797975, 0.14230026,
0.24348738, 0.33959405, 0.39992794, 0.77630536, 0.96993578])
In [75]:
len(s_sorted)
Out[75]:
10
Find the cumulative percentage (position) of each sorted values s_position = [0.1, 0.2, .....1] S_percentile = [-1.28, -0.84, ....]