editor_options: markdown: wrap: 72 —
This notebook walks through running and solving DEA models with
dea.py. You choose how many DMUs
(decision-making units), how many inputs/outputs, the
model (CCR or BCC) and the orientation
(input or output). Every code cell below was run and the outputs shown
are real.
Contents 1. Setup 2. The models in one page 3. Choose your parameters 4. Generate (or load) data 5. CCR input-oriented 6. CCR output-oriented 7. BCC input- and output-oriented 8. Slacks, peers and targets 9. Scale efficiency and returns to scale 10. A textbook check 11. Using your own CSV and the command line
Requirements: numpy, pandas,
scipy (the LPs are solved with
scipy.optimize.linprog using the HiGHS solver). Put
dea.py in your working folder.
pip install numpy pandas scipyfrom dea import generate_data, load_csv, solve_dea, scale_efficiency
import numpy as np, pandas as pdFor each DMU \(o\) with inputs \(x_o \in \mathbb{R}^m\) and outputs \(y_o \in \mathbb{R}^s\), and intensity weights \(\lambda \ge 0\) over all \(n\) DMUs:
Input-oriented (envelopment form)
\[\min_{\theta,\lambda}\ \theta \quad \text{s.t.}\quad \sum_j \lambda_j x_{ij} \le \theta\, x_{io},\quad \sum_j \lambda_j y_{rj} \ge y_{ro},\quad \lambda_j \ge 0\]
Efficiency \(=\theta^* \in (0,1]\): the proportion to which all inputs could be shrunk.
Output-oriented
\[\max_{\phi,\lambda}\ \phi \quad \text{s.t.}\quad \sum_j \lambda_j x_{ij} \le x_{io},\quad \sum_j \lambda_j y_{rj} \ge \phi\, y_{ro},\quad \lambda_j \ge 0\]
Efficiency \(=1/\phi^*\): \(\phi^*\ge 1\) is the factor by which all outputs could grow.
| Model | Returns to scale | Extra constraint |
|---|---|---|
| CCR (Charnes, Cooper & Rhodes 1978) | constant (CRS) | none |
| BCC (Banker, Charnes & Cooper 1984) | variable (VRS) | \(\sum_j \lambda_j = 1\) |
Phase 2. After \(\theta^*\) (or \(\phi^*\)) is found, a second LP fixes it and maximizes the sum of input slacks \(s^-\) and output slacks \(s^+\). A DMU is efficient (Pareto-Koopmans) only when the score is 1 and all slacks are 0. Score 1 with a positive slack is reported as weak.
Targets (projection onto the frontier): input-oriented \(\hat x = \theta^* x_o - s^-,\ \hat y = y_o + s^+\); output-oriented \(\hat x = x_o - s^-,\ \hat y = \phi^* y_o + s^+\).
Change these and re-run the notebook. MODEL is
"ccr" or "bcc"; ORIENTATION is
"input" or "output". A common rule of thumb is
\(n \ge \max\{m \cdot s,\ 3(m+s)\}\) so
the model has enough discrimination.
N_DMUS = 10
N_INPUTS = 2
N_OUTPUTS = 2
MODEL = "ccr" # "ccr" or "bcc"
ORIENTATION = "input" # "input" or "output"
SEED = 7
print("rule of thumb n >=", max(N_INPUTS * N_OUTPUTS, 3 * (N_INPUTS + N_OUTPUTS)))Output:
rule of thumb n >= 12generate_data creates strictly positive random data,
with outputs loosely tied to inputs so the result looks like a real
production process. Rows are DMUs.
X, Y, names = generate_data(N_DMUS, N_INPUTS, N_OUTPUTS, seed=SEED)
in_names = [f"x{i+1}" for i in range(N_INPUTS)]
out_names = [f"y{r+1}" for r in range(N_OUTPUTS)]
data = pd.DataFrame(np.hstack([X, Y]), index=names, columns=in_names + out_names)
print(data)Output:
x1 x2 y1 y2
DMU1 66.3 90.7 30.9 74.7
DMU2 79.8 30.3 33.7 43.2
DMU3 37.0 88.6 18.5 91.9
DMU4 10.5 83.9 25.1 96.5
DMU5 81.7 52.1 41.6 97.8
DMU6 37.3 35.1 19.9 39.0
DMU7 32.9 50.1 11.7 41.4
DMU8 55.4 59.8 37.8 58.2
DMU9 99.6 81.3 43.3 65.7
DMU10 66.0 99.0 60.4 77.9solve_dea solves one LP pair per DMU and returns a
DEAResult. summary() gives the efficiency
score, efficiency status, reference peers with their \(\lambda\) weights, and rank.
ccr_in = solve_dea(X, Y, model=MODEL, orientation=ORIENTATION,
names=names, input_names=in_names, output_names=out_names)
print(ccr_in.summary())Output:
efficiency theta ... peers rank
DMU ...
DMU1 0.6326 0.6326 ... DMU4(0.294), DMU5(0.363), DMU10(0.140) 9
DMU2 1.0000 1.0000 ... DMU2(1.000) 1
DMU3 0.7994 0.7994 ... DMU4(0.673), DMU5(0.276) 7
DMU4 1.0000 1.0000 ... DMU4(1.000) 1
DMU5 1.0000 1.0000 ... DMU5(1.000) 1
DMU6 0.8339 0.8339 ... DMU4(0.023), DMU5(0.269), DMU10(0.135) 6
DMU7 0.5828 0.5828 ... DMU4(0.220), DMU5(0.206) 10
DMU8 0.9460 0.9460 ... DMU2(0.077), DMU5(0.214), DMU10(0.435) 5
DMU9 0.7156 0.7156 ... DMU2(0.367), DMU5(0.225), DMU10(0.357) 8
DMU10 1.0000 1.0000 ... DMU10(1.000) 1
[10 rows x 5 columns]Reading it: a DMU with efficiency 0.80 could produce the same outputs with 80% of each input, by imitating the convex (here: non-negative) combination of its peers.
Under CRS the input- and output-oriented scores are reciprocals, so the efficiencies match the previous table; \(\phi\) tells you how much outputs could expand.
ccr_out = solve_dea(X, Y, "ccr", "output", names, in_names, out_names)
print(ccr_out.summary())
print("\nsame efficiency as input-oriented CCR:",
np.allclose(ccr_out.efficiency, solve_dea(X, Y, "ccr", "input").efficiency))Output:
efficiency phi ... peers rank
DMU ...
DMU1 0.6326 1.5809 ... DMU4(0.464), DMU5(0.573), DMU10(0.221) 9
DMU2 1.0000 1.0000 ... DMU2(1.000) 1
DMU3 0.7994 1.2509 ... DMU4(0.842), DMU5(0.345) 7
DMU4 1.0000 1.0000 ... DMU4(1.000) 1
DMU5 1.0000 1.0000 ... DMU5(1.000) 1
DMU6 0.8339 1.1992 ... DMU4(0.027), DMU5(0.322), DMU10(0.162) 6
DMU7 0.5828 1.7159 ... DMU4(0.377), DMU5(0.354) 10
DMU8 0.9460 1.0570 ... DMU2(0.082), DMU5(0.227), DMU10(0.460) 5
DMU9 0.7156 1.3973 ... DMU2(0.513), DMU5(0.315), DMU10(0.498) 8
DMU10 1.0000 1.0000 ... DMU10(1.000) 1
[10 rows x 5 columns]
same efficiency as input-oriented CCR: TrueBCC adds \(\sum\lambda_j = 1\), so each DMU is compared only with DMUs of a similar size. BCC scores are always \(\ge\) CCR scores, and input/output orientations can now differ.
bcc_in = solve_dea(X, Y, "bcc", "input", names, in_names, out_names)
bcc_out = solve_dea(X, Y, "bcc", "output", names, in_names, out_names)
compare = pd.DataFrame({
"CCR": ccr_in.efficiency.round(4),
"BCC-in": bcc_in.efficiency.round(4),
"BCC-out": bcc_out.efficiency.round(4),
}, index=names)
print(compare)Output:
CCR BCC-in BCC-out
DMU1 0.6326 0.6596 0.7778
DMU2 1.0000 1.0000 1.0000
DMU3 0.7994 0.8128 0.9476
DMU4 1.0000 1.0000 1.0000
DMU5 1.0000 1.0000 1.0000
DMU6 0.8339 1.0000 1.0000
DMU7 0.5828 0.9365 0.6965
DMU8 0.9460 0.9964 0.9955
DMU9 0.7156 0.7217 0.8114
DMU10 1.0000 1.0000 1.0000Slacks are the extra, non-radial improvements left after the proportional reduction. Targets are where each DMU lands on the frontier.
print(bcc_in.slacks())Output:
s-_x1 s-_x2 s+_y1 s+_y2
DMU
DMU1 0.0 0.0 0.0000 0.0000
DMU2 0.0 0.0 0.0000 0.0000
DMU3 0.0 0.0 10.1602 0.0000
DMU4 0.0 0.0 0.0000 0.0000
DMU5 0.0 0.0 0.0000 0.0000
DMU6 0.0 0.0 0.0000 0.0000
DMU7 0.0 0.0 9.4592 11.5243
DMU8 0.0 0.0 0.0000 0.0000
DMU9 0.0 0.0 0.0000 0.0000
DMU10 0.0 0.0 0.0000 0.0000print(bcc_in.targets())Output:
x1* x2* y1* y2*
DMU
DMU1 43.728 59.822 30.900 74.700
DMU2 79.800 30.300 33.700 43.200
DMU3 30.072 72.010 28.660 91.900
DMU4 10.500 83.900 25.100 96.500
DMU5 81.700 52.100 41.600 97.800
DMU6 37.300 35.100 19.900 39.000
DMU7 30.810 46.917 21.159 52.924
DMU8 55.199 59.583 37.800 58.200
DMU9 71.879 58.672 43.300 65.700
DMU10 66.000 99.000 60.400 77.900The full \(\lambda\) matrix (rows = evaluated DMU, columns = peer) is available too:
print(bcc_in.lambdas_df())Output:
DMU1 DMU2 DMU3 DMU4 DMU5 DMU6 DMU7 DMU8 DMU9 DMU10
DMU
DMU1 0.0 0.0000 0.0 0.2914 0.2590 0.3541 0.0 0.0 0.0 0.0954
DMU2 0.0 1.0000 0.0 0.0000 0.0000 0.0000 0.0 0.0 0.0 0.0000
DMU3 0.0 0.0000 0.0 0.6718 0.2427 0.0855 0.0 0.0 0.0 0.0000
DMU4 0.0 0.0000 0.0 1.0000 0.0000 0.0000 0.0 0.0 0.0 0.0000
DMU5 0.0 0.0000 0.0 0.0000 1.0000 0.0000 0.0 0.0 0.0 0.0000
DMU6 0.0 0.0000 0.0 0.0000 0.0000 1.0000 0.0 0.0 0.0 0.0000
DMU7 0.0 0.0000 0.0 0.2422 0.0000 0.7578 0.0 0.0 0.0 0.0000
DMU8 0.0 0.0928 0.0 0.0000 0.0750 0.4620 0.0 0.0 0.0 0.3702
DMU9 0.0 0.3725 0.0 0.0000 0.2001 0.0838 0.0 0.0 0.0 0.3436
DMU10 0.0 0.0000 0.0 0.0000 0.0000 0.0000 0.0 0.0 0.0 1.0000# peers of one DMU as a dict
worst = int(np.argmin(bcc_in.efficiency))
print(names[worst], "->", bcc_in.peers(worst))Output:
DMU1 -> {'DMU4': 0.2914, 'DMU5': 0.259, 'DMU6': 0.3541, 'DMU10': 0.0954}Technical efficiency (CCR) splits into pure technical efficiency (BCC) times scale efficiency: \(SE = \text{CCR} / \text{BCC}\). Returns to scale are read from \(\sum\lambda_j\) in the CCR solution: \(<1\) increasing (IRS), \(=1\) constant (CRS), \(>1\) decreasing (DRS). The CCR \(\lambda\) may not be unique, so treat the RTS label as indicative for DMUs with alternative optima.
print(scale_efficiency(X, Y, "input", names, in_names, out_names))Output:
CCR (TE) BCC (PTE) SE sum_lambda RTS
DMU
DMU1 0.6326 0.6596 0.9591 0.7961 IRS
DMU2 1.0000 1.0000 1.0000 1.0000 CRS
DMU3 0.7994 0.8128 0.9836 0.9486 IRS
DMU4 1.0000 1.0000 1.0000 1.0000 CRS
DMU5 1.0000 1.0000 1.0000 1.0000 CRS
DMU6 0.8339 1.0000 0.8339 0.4265 IRS
DMU7 0.5828 0.9365 0.6223 0.4262 IRS
DMU8 0.9460 0.9964 0.9495 0.7269 IRS
DMU9 0.7156 0.7217 0.9916 0.9495 IRS
DMU10 1.0000 1.0000 1.0000 1.0000 CRSCooper, Seiford & Tone (Data Envelopment Analysis, ch. 1) use eight stores with one input (employees) and one output (sales). Their CCR scores are A 0.5, B 1, C 0.667, D 0.75, E 0.8, F 0.4, G 0.5, H 0.625, and BCC input-oriented C 0.833, F 0.5, with A, B, E, H on the VRS frontier. The solver reproduces them:
Xs = [[2],[3],[3],[4],[5],[5],[6],[8]]
Ys = [[1],[3],[2],[3],[4],[2],[3],[5]]
stores = list("ABCDEFGH")
print(pd.DataFrame({
"CCR": solve_dea(Xs, Ys, "ccr", "input", stores).efficiency.round(4),
"BCC-in": solve_dea(Xs, Ys, "bcc", "input", stores).efficiency.round(4),
"BCC-out":solve_dea(Xs, Ys, "bcc", "output", stores).efficiency.round(4),
}, index=stores))Output:
CCR BCC-in BCC-out
A 0.5000 1.0000 1.0000
B 1.0000 1.0000 1.0000
C 0.6667 0.8333 0.6667
D 0.7500 0.7500 0.8571
E 0.8000 1.0000 1.0000
F 0.4000 0.5000 0.5000
G 0.5000 0.5000 0.6923
H 0.6250 1.0000 1.0000One row per DMU; pick input and output columns by name.
csv = "example_stores.csv"
pd.DataFrame({"Store": list("ABCDEFGH"),
"Employees": [4, 7, 8, 4, 2, 10, 3, 5],
"FloorArea": [3, 3, 1, 2, 4, 1, 7, 5],
"Sales": [1, 1, 1, 1, 1, 1, 1, 1]}).to_csv(csv, index=False)
X2, Y2, n2 = load_csv(csv, inputs=["Employees", "FloorArea"], outputs=["Sales"], name_col="Store")
res = solve_dea(X2, Y2, "ccr", "input", n2, ["Employees", "FloorArea"], ["Sales"])
print(res.summary())
print()
print(res.slacks())Output:
efficiency theta efficient peers rank
DMU
A 0.8571 0.8571 no D(0.714), E(0.286) 5
B 0.6316 0.6316 no C(0.105), D(0.895) 7
C 1.0000 1.0000 yes C(1.000) 1
D 1.0000 1.0000 yes D(1.000) 1
E 1.0000 1.0000 yes E(1.000) 1
F 1.0000 1.0000 weak C(1.000) 1
G 0.6667 0.6667 no E(1.000) 6
H 0.6000 0.6000 no D(0.500), E(0.500) 8
s-_Employees s-_FloorArea s+_Sales
DMU
A 0.0 0.0000 0.0
B 0.0 0.0000 0.0
C 0.0 0.0000 0.0
D 0.0 0.0000 0.0
E 0.0 0.0000 0.0
F 2.0 0.0000 0.0
G 0.0 0.6667 0.0
H 0.0 0.0000 0.0Store F scores 1 but carries 2 surplus employees, so it is only weakly efficient; G has a floor-area slack on top of its radial inefficiency.
Command-line equivalents:
# you choose the number of DMUs and the model
python dea.py --n-dmus 15 --n-inputs 3 --n-outputs 2 --model bcc --orientation output
# CCR + BCC + scale efficiency in one go, saved to CSV
python dea.py --n-dmus 20 --model both --save results.csv
# your own data
python dea.py --data example_stores.csv --inputs Employees FloorArea --outputs Sales \\
--name-col Store --model ccr --orientation inputNotes and limits - CCR/BCC need strictly positive inputs and outputs; the solver raises an error otherwise. - Performance: 200 DMUs × 5 variables solves in about a second. - Ties in efficiency share the same rank; efficient DMUs all rank 1 (use super-efficiency if you need to break those ties).