Numerical and Statistical Programming

Introduction

Numerical and statistical programming provides the computational foundation for business analytics. It enables analysts to transform data into mathematical objects, perform calculations at scale, quantify uncertainty, test claims, and build predictive models. In Python, this work is commonly conducted with NumPy for numerical arrays, pandas for labelled tabular data, SciPy for scientific and statistical procedures, and statsmodels or scikit-learn for statistical modelling and machine learning.

The central principle is vectorized computation: express an operation on an entire array rather than writing an explicit loop for every observation. Vectorization improves clarity, reduces programming errors, and usually exploits optimized low-level implementations. However, computational speed does not guarantee statistical validity. A statistically inappropriate model can produce precise but misleading results. Professional analytics therefore combines programming discipline, mathematical reasoning, assumption checking, and responsible interpretation.

Learning Objectives

After completing this chapter, a learner should be able to:

1. Represent numerical data using vectors, matrices, and data frames.
2. Apply vectorized operations, broadcasting, masking, and aggregation.
3. Perform matrix operations relevant to analytics and regression.
4. Simulate random variables and use simulation to study probability.
5. Calculate and interpret descriptive statistics and confidence intervals.
6. Design samples and conduct hypothesis tests responsibly.
7. Fit, diagnose, and interpret regression models.
8. Recognize computational, statistical, ethical, and reproducibility risks.

1. Numerical Objects and Array Programming

1.1 Scalars, vectors, matrices, and data frames

A scalar is a single number. A vector is an ordered one-dimensional collection, and a matrix is a rectangular two-dimensional collection. An array generalizes matrices to any number of dimensions. A data frame is a labelled table whose columns may have different data types.

In NumPy, an array has a shape, data type, and typically homogeneous storage. For example:

import numpy as np

sales = np.array([120, 150, 135, 170])
prices = np.array([10.0, 9.5, 11.0, 10.5])
revenue = sales * prices

The expression sales * prices performs elementwise multiplication, producing [1200, 1425, 1485, 1785]. It is not matrix multiplication. Matrix multiplication uses @ or np.matmul.

A vector x with n observations can be written as:

The superscript T denotes transpose. Shape is operationally important: a vector with shape (n,), a column matrix with shape (n, 1), and a row matrix with shape (1, n) may behave differently in broadcasting and matrix multiplication.

1.2 Indexing, slicing, and masks

Python uses zero-based indexing. Slicing selects ranges without copying in every context, while Boolean masks select observations satisfying a condition:

high_value = revenue[revenue > 1400]

Masks are fundamental for data quality rules, segment analysis, and outlier review. A mask should be explicit about missing values and business meaning. For example, revenue > 1400 does not mean that an observation is anomalous; it only applies a threshold.

1.3 Vectorization and broadcasting

A vectorized transformation applies one operation to many values:

z = (x - x.mean()) / x.std(ddof=1)

Broadcasting permits arrays of compatible shapes to interact. If a matrix has shape (m, p) and a vector of column means has shape (p,), subtracting the vector from the matrix subtracts each column mean from the corresponding column. Broadcasting avoids repeated code but requires careful shape checking.

Useful operations include sum, mean, median, min, max, std, argmax, and where. The axis argument determines the direction of aggregation. For a matrix with rows as observations and columns as variables, axis=0 usually produces one result per column, whereas axis=1 produces one result per row.

1.4 Missing values and numerical precision

Missing values must not be silently treated as zero. NumPy commonly represents missing floating-point values as np.nan; pandas provides methods such as isna, fillna, and dropna. The correct response depends on the data-generating process and the analysis objective.

Floating-point numbers are approximations. Thus, direct equality tests such as a == b can fail after arithmetic. Use a tolerance where appropriate:

np.isclose(a, b, rtol=1e-9, atol=1e-12)

Large arrays also raise memory and performance concerns. Appropriate dtypes, chunking, sparse matrices, and avoidance of unnecessary copies can materially improve performance.

2. Matrix Algebra for Analytics

Matrices provide a compact language for multivariate data, optimization, and regression. Let A be an m × n matrix and B an n × p matrix. Their product AB is an m × p matrix, with elements:

The inner dimensions must agree. Matrix multiplication is not commutative: AB generally differs from BA.

The dot product of vectors a and b of equal length is:

The Euclidean norm is:

For a square matrix A, an inverse A⁻¹ exists only when A is nonsingular. Although the formal solution to Ax = b is x = A⁻¹b, computing an explicit inverse is often less stable and slower than solving the system directly:

x = np.linalg.solve(A, b)

Important diagnostics include the determinant, rank, eigenvalues, and condition number. A high condition number indicates that small input errors may create large changes in the solution. In regression, multicollinearity can make coefficient estimates unstable even when predictions remain adequate.

2.1 Covariance and correlation matrices

For observations of p variables, the sample covariance between variables X and Y is:

The sample standard deviation is:

The sample correlation is:

These formulas assume paired observations and finite, appropriately handled data. Covariance has the product of the variables’ units; correlation is dimensionless and lies between −1 and 1. Correlation measures linear association, not causation.

3. Probability and Simulation

Probability programming is useful when analytical calculations are difficult, when uncertainty must be communicated, or when a proposed decision rule needs stress testing. A random variable X has a probability distribution describing possible values and their probabilities.

The expected value and variance of a discrete random variable are:

where μ = E[X]. For a continuous variable, sums are replaced by integrals over the density, provided the relevant expectation exists.

A Bernoulli variable takes value 1 with probability p and 0 with probability 1−p. A binomial variable counts successes in n independent Bernoulli trials with common success probability p:

for k = 0, 1, …, n. The assumptions of independence and constant p are substantive, not merely programming details.

3.1 Monte Carlo simulation

Monte Carlo simulation repeatedly generates random inputs and records outputs. A reproducible workflow uses a local random-number generator:

rng = np.random.default_rng(42)
arrivals = rng.poisson(lam=4, size=100_000)
probability = np.mean(arrivals >= 6)

The simulation estimate of a probability is a sample proportion. If the true probability is q and N independent simulations are used, its standard error is:

Since q is unknown, q̂ may be substituted for q. Increasing N generally reduces Monte Carlo error at a rate proportional to 1/√N; it does not correct a misspecified distribution.

Random seeds make results reproducible, but a seed is not a substitute for documenting distributions, parameters, code versions, and data. In production, random-number streams should be managed so parallel processes do not unintentionally duplicate draws.

4. Descriptive Statistical Analysis

Descriptive statistics summarize observed data; they do not automatically generalize beyond it. Measures of central tendency include the mean, median, and mode. The mean is:

It is sensitive to extreme values. The median is the middle ordered observation, or the average of the two middle observations for an even sample. Measures of spread include range, interquartile range, variance, standard deviation, and median absolute deviation.

The sample variance uses n−1 in the denominator:

The n−1 correction is appropriate when estimating a population variance from an independently sampled population under the usual finite-variance assumptions. It is not universally required: if describing the complete finite dataset, division by n may be the descriptive convention.

Quantiles divide ordered data according to cumulative proportion. The interquartile range is Q₃−Q₁. Box plots, histograms, empirical cumulative distributions, scatter plots, and grouped summaries help reveal skewness, multimodality, dependence, and data errors. Visualizations should show units, denominators, missing-data treatment, and meaningful labels.

A standardized value is:

It indicates distance from the sample mean in sample-standard-deviation units. Standardization does not make a variable normally distributed.

5. Sampling and Statistical Inference

A population is the full set of units of interest; a sample is the observed subset. Probability sampling methods include simple random, stratified, cluster, and systematic sampling. Convenience sampling can be useful for exploratory work but generally cannot support unbiased population inference without strong assumptions.

A sampling frame should cover the target population adequately. Undercoverage, nonresponse, voluntary participation, and measurement error can create bias that a larger sample does not eliminate. Stratification can improve precision when strata are internally similar and meaningfully represented.

5.1 Sampling distributions and the central limit theorem

A statistic is a function of sample data. Repeated samples produce a sampling distribution. Under suitable conditions—often independence, finite variance, and a sufficiently large sample—the central limit theorem implies that the sample mean is approximately normal:

Here μ and σ are the population mean and standard deviation, n is sample size, and N(0,1) is the standard normal distribution. The approximation can be poor for small samples or highly skewed populations.

If σ is unknown, a confidence interval for a mean commonly uses the t distribution:

The interval assumes a random sample and, for small samples, approximate normality of the population or robustness supported by the design. A 95% confidence procedure means that, over repeated samples and under its assumptions, approximately 95% of constructed intervals contain the fixed population mean. It does not mean there is a 95% probability that one fixed parameter lies in this particular interval under frequentist interpretation.

5.2 Hypothesis testing

A hypothesis test begins with a null hypothesis H₀ and an alternative H₁, a test statistic, a decision rule, and an error criterion. The p-value is the probability, assuming H₀ and the model are true, of obtaining a result at least as extreme as the observed result in the specified direction.

A small p-value provides evidence against H₀; it does not measure the probability that H₀ is true, the size of an effect, or its practical importance. Type I error rejects a true H₀; Type II error fails to reject a false H₀. Statistical power is 1−β, where β is the Type II error probability under a specified alternative.

For a one-sample mean test with known σ:

where μ₀ is the null value. If σ is estimated, a one-sample t statistic is used:

The test choice must match the design. Paired observations require a paired analysis of within-pair differences; independent groups require an independent-samples method. Multiple testing inflates the chance of false positives. Corrections such as Bonferroni control a family-wise error rate but may reduce power; false-discovery-rate procedures address a different objective.

Effect sizes and confidence intervals should accompany p-values. For two independent means, a standardized effect may be expressed as Cohen’s d, with the exact denominator choice documented. Practical significance depends on costs, benefits, and domain thresholds, not solely on a conventional α such as 0.05.

6. Regression Programming

Regression models describe conditional relationships and support prediction or explanation. A multiple linear regression is:

Here yᵢ is the response, xᵢⱼ are predictors, β₀ is the intercept, βⱼ are partial slope coefficients, and εᵢ is an error term. In matrix form:

The ordinary least-squares estimator minimizes the sum of squared residuals:

If X has full column rank, the algebraic solution is:

In numerical software, QR or singular-value-decomposition methods are preferred to explicitly calculating the inverse. The classical inferential interpretation commonly assumes linearity in parameters, independent errors, exogeneity, constant error variance for conventional standard errors, and—especially in small samples—approximately normal errors. These assumptions concern the model and design, not merely the dependent variable.

The fitted value is ŷᵢ and the residual is eᵢ = yᵢ−ŷᵢ. The coefficient of determination is:

R² is the proportion of sample variation explained relative to an intercept-only model. It is not a universal measure of predictive quality and can increase when irrelevant predictors are added. Adjusted R² penalizes model size, but cross-validation and out-of-sample metrics are often more informative for prediction.

A typical implementation is:

import statsmodels.api as sm

X = df[["advertising", "price"]]
X = sm.add_constant(X)
y = df["sales"]
model = sm.OLS(y, X, missing="drop").fit()
print(model.summary())

Diagnostics include residual-versus-fitted plots for nonlinearity and unequal variance, Q–Q plots for distributional assessment, leverage and influence measures for unusual observations, and variance-inflation diagnostics for collinearity. Robust standard errors can address some forms of heteroscedasticity, but they do not repair omitted-variable bias, reverse causality, poor measurement, or a nonrepresentative sample.

Categorical variables require a reference category and appropriate encoding. Interaction terms represent changing effects; for example, an interaction between price and customer segment means the price slope differs by segment. Extrapolation beyond the observed predictor range is especially risky.

7. Applied Case: Simulating and Evaluating a Service Campaign

The following is an illustrative case.

Background and problem

A subscription firm wants to evaluate whether a reminder campaign increases monthly renewal. Management observes that customers receiving reminders renew more often, but reminder assignment was historically targeted toward customers already likely to renew. The observed association may therefore reflect selection bias.

Analysis and approach

The analytics team first summarizes renewal rates by campaign status, checks missingness, and compares customer tenure and prior engagement. It then creates a reproducible simulation of a hypothetical randomized experiment: for each of 100,000 simulated customers, baseline renewal probability is drawn from a specified logistic model based on tenure and engagement; treatment assignment is generated independently with probability 0.5; treatment adds a pre-specified log-odds effect.

For each simulated experiment, the team calculates the difference in renewal proportions and records whether a 95% confidence interval excludes zero. Repeating the experiment estimates sampling variability and statistical power under the assumed effect. The team also fits a logistic regression to the simulated outcomes, including treatment and baseline covariates.

The actual evaluation then uses a randomized rollout where feasible. If randomization is impossible, analysts use a clearly stated observational design, such as covariate adjustment or matching, while emphasizing that conditional association is not automatically causal.

Outcome and interpretation

Simulation shows that a modest campaign effect can be missed in small samples, especially when baseline renewal is variable. The randomized rollout provides an unbiased treatment comparison in expectation, subject to implementation compliance and attrition. Management receives the estimated effect, confidence interval, incremental renewals, campaign cost, and sensitivity analysis rather than a binary “significant/not significant” label.

Learning points

1. Simulation quantifies uncertainty but depends on its assumed data-generating process.
2. Random assignment addresses confounding more directly than regression adjustment alone.
3. Statistical significance must be combined with effect size and economic value.
4. Reproducible seeds, documented assumptions, and versioned code are part of professional analysis.

8. Best Practices, Limitations, and Governance

Use vectorized operations but inspect shapes and units. Separate data preparation, analysis, visualization, and reporting. Validate inputs, preserve raw data, record transformations, and test functions with known examples. Use stable numerical solvers, avoid unnecessary matrix inversion, and monitor overflow, underflow, and loss of precision.

Statistical best practice requires defining the estimand before analysis, documenting inclusion criteria, distinguishing exploratory from confirmatory work, and reporting uncertainty. Do not select a test after inspecting results, repeatedly test many outcomes without disclosure, or present a predictive model as causal evidence.

Privacy and governance are also relevant. Financial, health, employee, and customer data may be regulated or contractually restricted. Analysts should apply data minimization, access controls, secure storage, de-identification where appropriate, and retention policies. Reproducibility must not mean exposing confidential records. Synthetic data can support development, but it does not automatically preserve the statistical or privacy properties of real data.

Contemporary workflows increasingly use automatic differentiation, probabilistic programming, distributed arrays, and AI-assisted coding. These tools can accelerate analysis, but generated code and generated interpretations require human review, tests, provenance, and disclosure. Libraries should be pinned or recorded, random seeds managed, and environments reproducible. A notebook is not by itself a reproducible study unless dependencies, inputs, execution order, and outputs are controlled.

Chapter Summary

Numerical and statistical programming links mathematical ideas with reliable analytical practice. Arrays, vectorization, broadcasting, masks, and aggregation provide efficient computational building blocks. Matrix algebra supports multivariate analysis and regression, while careful numerical methods avoid instability. Simulation translates probability assumptions into empirical evidence and helps quantify Monte Carlo error. Descriptive statistics summarize data, whereas sampling theory, confidence intervals, and hypothesis tests support cautious generalization. Regression estimates conditional relationships and requires diagnostic assessment of design and assumptions. Across all methods, interpretation depends on data quality, sampling, model specification, uncertainty, and ethical governance. Professional analysts therefore report assumptions, effect sizes, uncertainty, limitations, and reproducible computational procedures.

Key Terms

Vectorization: Applying an operation to an entire array rather than iterating observation by observation.

Broadcasting: Automatic alignment of compatible array shapes for elementwise operations.

Mask: A Boolean array used to select observations.

Matrix multiplication: Product of matrices satisfying compatible inner dimensions.

Condition number: A measure of sensitivity of a numerical problem to input perturbations.

Covariance: Joint measure of variation between two variables.

Correlation: Standardized measure of linear association.

Monte Carlo simulation: Repeated random sampling to approximate a probability or quantity.

Sampling frame: Operational list or mechanism identifying units available for sampling.

Sampling distribution: Distribution of a statistic over repeated samples.

Central limit theorem: Result describing approximate normality of suitably standardized sums or means under conditions.

Confidence interval: Interval produced by a procedure designed to cover a parameter at a stated long-run rate.

Null hypothesis: Default claim evaluated by a statistical test.

p-value: Probability, under the null model, of a result at least as extreme as observed.

Type I error: Rejection of a true null hypothesis.

Statistical power: Probability of rejecting a false null under a specified alternative.

Residual: Observed response minus fitted response.

Heteroscedasticity: Nonconstant variance of model errors.

Multicollinearity: Strong dependence among predictors that can destabilize coefficient estimates.

Reproducibility: Ability to obtain consistent analytical results using documented data, code, and environment.

Review Questions

Short Answer Questions

1. Distinguish elementwise multiplication from matrix multiplication in NumPy.
2. What does a high matrix condition number imply?
3. Why is a random seed useful in simulation?
4. Explain the difference between a p-value and an effect size.
5. State two assumptions relevant to ordinary least-squares inference.

Long Answer Questions

1. Explain vectorization, broadcasting, masking, and aggregation using a business dataset.
2. Derive and interpret the sample covariance and correlation formulas, including their units and limitations.
3. Describe how a confidence interval should be interpreted and explain common misinterpretations.
4. Explain the OLS regression objective, matrix solution, and reasons to avoid explicit matrix inversion.
5. Discuss how sampling bias, multiple testing, and model diagnostics affect managerial conclusions.

Analytical / Application Questions

1. A retailer simulates 50,000 days of demand from a Poisson distribution but obtains a service-level estimate that changes substantially across runs. Design a reproducibility and Monte Carlo-error assessment.
2. An advertising coefficient is positive and statistically significant in an observational regression. Identify at least four reasons why this result may not establish a causal advertising effect, and propose a stronger evaluation design.
3. A customer dataset contains missing income values, extreme purchase amounts, and highly correlated predictors. Develop a defensible workflow covering data treatment, descriptive analysis, regression diagnostics, and reporting.

Multiple Choice Questions with Answers

Question 1

Which NumPy expression performs elementwise multiplication of two compatible arrays x and y?

A. x * y

B. x @ y

C. np.linalg.solve(x, y)

D. x.T + y

Correct Answer: A. x * y

Question 2

For a matrix A with shape (m, n) and B with shape (n, p), what is the shape of A @ B?

A. (m, p)

B. (n, n)

C. (m, n)

D. (p, m)

Correct Answer: A. (m, p)

Question 3

What does a high condition number generally indicate?

A. The matrix contains only zeros

B. The numerical solution may be highly sensitive to input errors

C. The matrix must be symmetric

D. The sample size is large

Correct Answer: B. The numerical solution may be highly sensitive to input errors

Question 4

If 10,000 independent simulations estimate a probability as 0.20, the approximate plug-in Monte Carlo standard error is closest to:

A. 0.0004

B. 0.0040

C. 0.0126

D. 0.1265

Correct Answer: C. 0.0126

Question 5

Which statement about a 95% frequentist confidence interval is correct?

A. There is a 95% probability that the fixed parameter lies in this computed interval

B. It contains 95% of the observations

C. The procedure captures the parameter in approximately 95% of repeated samples under its assumptions

D. It proves the null hypothesis is false

Correct Answer: C. The procedure captures the parameter in approximately 95% of repeated samples under its assumptions

Question 6

A p-value is best described as:

A. The probability that the null hypothesis is true

B. The probability, under the null model, of a result at least as extreme as observed

C. The probability that the alternative hypothesis is true

D. The practical size of an effect

Correct Answer: B. The probability, under the null model, of a result at least as extreme as observed

Question 7

In ordinary least squares, the fitted coefficients minimize:

A. The sum of absolute predictor values

B. The sum of squared residuals

C. The number of observations

D. The sample correlation only

Correct Answer: B. The sum of squared residuals

Question 8

Which diagnostic is especially useful for identifying nonconstant residual variance?

A. Residual-versus-fitted plot

B. A random-number seed

C. A covariance matrix alone

D. An index lookup

Correct Answer: A. Residual-versus-fitted plot

Question 9

If a sample mean is 80, the null mean is 75, the sample standard deviation is 10, and n = 25, the one-sample t statistic is:

A. 0.50

B. 1.00

C. 2.50

D. 12.50

Correct Answer: C. 2.50

Question 10

Which design most directly supports a causal comparison between a campaign and renewal?

A. Comparing volunteers with nonvolunteers

B. Comparing historical averages only

C. Randomly assigning eligible customers to campaign and control groups

D. Selecting only customers with high prior renewal rates

Correct Answer: C. Randomly assigning eligible customers to campaign and control groups