5 Scipy.stats-Methods zur Simulation von „Was-wäre-wenn“-Szenarien

# Einführung

Daten sind selten statisch. Entscheidungen sind selten risikofrei. Als Datenwissenschaftler werden Sie häufig gebeten, Geschäftsannahmen einem Stresstest zu unterziehen, Verteilungsunsicherheiten zu untersuchen oder various Realitäten zu simulieren.

• „Was wäre, wenn unsere tägliche aktive Nutzerakquise das Doppelte kosten würde?“
• „Was passiert, wenn unser Serververkehr während einer Werbeveranstaltung um 300 % ansteigt?“
• „Wie hoch ist die Wahrscheinlichkeit, dass unsere Betriebsverluste in diesem Quartal 50.000 US-Greenback übersteigen?“

Beantwortung dieser Fragen was ist, wenn Fragen erfordert den Übergang von einfachen Punktschätzungen (wie dem einfachen Mittelwert) zu robustem, probabilistischem Denken. Während viele Praktiker sofort auf umfangreiche Simulations-Engines umsteigen, enthält der standardmäßige wissenschaftliche Python-Stack bereits ein nicht ausreichend genutztes Arbeitstier für genau diese Artwork der Modellierung: scipy.stats. Abgesehen von seinem allgemeinen Ruf für die Durchführung einfacher Hypothesentests oder die Berechnung von p-Werten, scipy.stats Bietet eine einheitliche programmatische Schnittstelle zum Parametrisieren, Abtasten und Berechnen von Risikometriken für Dutzende kontinuierlicher und diskreter Wahrscheinlichkeitsverteilungen.

In diesem Artikel werfen wir einen Blick unter die Haube scipy.statsDabei werden fünf wesentliche Methods untersucht, um leistungsstarke, strenge Simulationen nur mit NumPy und SciPy zu entwerfen.

# 1. Einfrieren von Verteilungen zur Parametrisierung von Szenarien

Bei der Modellierung von Szenarien möchten Sie häufig unterschiedliche Zustände der Welt darstellen: eine konservative Basislinie, einen optimistischen Finest-Case und einen pessimistischen Worst-Case. Im Standardprozedurcode könnten Sie diese darstellen, indem Sie Wörterbücher mit Parametern (z. B. Standort) mit sich herumtragen loc und Maßstab scale) und entpacken Sie sie jedes Mal in Funktionen, wenn Sie eine Wahrscheinlichkeit bewerten oder eine Stichprobe ziehen müssen.

Ein überlegenes, objektorientiertes Muster ist Einfrieren Verteilungen. In scipy.statsAufruf einer Verteilungsklasse (z. B stats.norm, stats.lognormoder stats.gamma) und die direkte Übergabe Ihrer Parameter an den Konstruktor gibt eine „eingefrorene“ Zufallsvariable zurück (eine Instanz von rv_frozen).

Dieses eingeschlossene Objekt kapselt das gesamte Wahrscheinlichkeitsmodell. Sie können es als einzelnes Objekt in Ihrer Codebasis weitergeben, Szenariodefinitionen im Handumdrehen austauschen und Methoden wie aufrufen .rvs(), .pdf(), .cdf()oder .ppf() ohne jemals die Parameter erneut angeben zu müssen.

Angenommen, wir modellieren die tägliche Produktnachfrage. Bei einer manuellen Implementierung müssen wir Wörterbücher oder Variablen durch jede Section unserer Verarbeitungspipeline ziehen:

import scipy.stats as stats

# Defining uncooked situation parameters
situations = {
    "recession": {"imply": 800, "std": 250},
    "baseline": {"imply": 1200, "std": 150},
    "growth": {"imply": 1800, "std": 300}
}

# Drawing samples or evaluating threat requires handbook parameter unpacking
def simulate_demand(scenario_name, measurement=5):
    params = situations(scenario_name)
    # Passing parameters inside each statistical name
    samples = stats.norm.rvs(loc=params("imply"), scale=params("std"), measurement=measurement)
    p_exceed_1000 = 1 - stats.norm.cdf(1000, loc=params("imply"), scale=params("std"))
    return samples, p_exceed_1000

for title in situations.keys():
    samples, prob = simulate_demand(title)
    print(f"{title.capitalize()} -> Prob > 1000: {prob:.2%}")

Hier ist eine elegantere Various. Durch die Instanziierung der Verteilung friert SciPy die Parameter ein und gibt uns ein in sich geschlossenes, sauberes Szenarioobjekt:

import scipy.stats as stats

# With scipy, freeze the distribution parameters right into a named object
scenario_models = {
    "recession": stats.norm(loc=800, scale=250),
    "baseline": stats.norm(loc=1200, scale=150),
    "growth": stats.norm(loc=1800, scale=300)
}

# The pipeline merely accepts a frozen RV object, decoupling logic from parameters
def run_scenario_pipeline(rv_frozen, measurement=5):

    # Lock-in strategies are known as instantly on the item
    samples = rv_frozen.rvs(measurement=measurement)

    # sf() is survival perform (1 - CDF)
    p_exceed_1000 = rv_frozen.sf(1000)

    return samples, p_exceed_1000

for title, mannequin in scenario_models.objects():
    _, prob = run_scenario_pipeline(mannequin)
    print(f"{title.capitalize()} State of affairs | Prob demand > 1000: {prob:.2%}")

Ausgabe:

Recession State of affairs | Prob demand > 1000: 21.19%
Baseline State of affairs | Prob demand > 1000: 90.88%
Increase State of affairs | Prob demand > 1000: 99.62%

Indem wir unsere Parameter einfrieren, trennen wir unsere mathematischen Annahmen von unserer Ausführungslogik. Wenn wir unser Nachfragemodell von einer Normalverteilung auf eine verzerrte logarithmische Normalverteilung umstellen möchten, müssen wir nur die einzelne Zeile ändern, in der wir das eingefrorene Objekt initialisieren. Die nachgelagerte Pipeline bleibt unberührt.

# 2. Monte-Carlo-Simulation mit .rvs() zur Unsicherheitsquantifizierung

Bei der Geschäftsplanung erstellen Praktiker häufig Tabellenkalkulationen, die den Umsatz anhand statischer Punktschätzungen berechnen – z income = every day site visitors * conversion charge * common order worth. Allerdings ist jede dieser Eingaben höchst unsicher. Durch die Multiplikation von Durchschnittswerten wird die zusammengesetzte Varianz ausgeblendet, was zur Folge hat Fehler der Durchschnittswerteoder die systematische Unterschätzung des Risikos.

Um diese Unsicherheit zu quantifizieren, können wir verwenden Monte-Carlo-Simulation. Anstatt statische Zahlen anzunehmen, stellen wir jede Variable als Wahrscheinlichkeitsverteilung dar. Mit der .rvs(measurement=n) Mit der Methode für unsere eingefrorenen Verteilungen können wir sofort $N$ Zufallsstichproben aus allen Eingaben ziehen, sie durch unsere eigene Formel in einer vektorisierten NumPy-Pipeline weitergeben und die vollständige Wahrscheinlichkeitsverteilung des Endergebnisses auswerten.

Bei der Berechnung eines statischen besten/schlechtesten Falls werden die gemeinsame Wahrscheinlichkeit von Ereignissen und die tatsächliche Verteilung der Ergebnisse nicht berücksichtigt, während das Schreiben manueller Schleifen in reinem Python unglaublich langsam ist.

import random

# Brittle level estimates
avg_traffic = 50000
avg_conversion = 0.05
avg_aov = 60.0

expected_revenue = avg_traffic * avg_conversion * avg_aov
print(f"Level Estimate Anticipated Income: ${expected_revenue:,.2f}")

# Sluggish, iterative handbook sampling loop
n_sims = 100000
revenue_clunky = ()
for _ in vary(n_sims):
    # Simulating regular site visitors and uniform conversion charges
    site visitors = random.gauss(50000, 5000)
    conversion = random.uniform(0.04, 0.06)
    aov = random.gammavariate(15, 4.0)
    revenue_clunky.append(site visitors * conversion * aov)

Ausgabe:

Level Estimate Anticipated Income: $150,000.00

Durch die Nutzung scipy.stats Verteilungsmodelle und gleichzeitige Stichprobenziehung mit .rvs()können wir die gesamte Simulation in einer einzigen vektorisierten NumPy-Operation berechnen. Lassen Sie uns das Downside in vier verschiedenen Schritten angehen:

1. Definieren Sie geeignete Verteilungsformen für unser Szenario
2. Zeichnen Sie N Proben
3. Ausbreitung durch Geschäftslogik (über Vektorisierung)
4. Extrahieren Sie Risikoperzentile

import numpy as np
import scipy.stats as stats

# 1. Outline applicable distribution shapes for our situation
np.random.seed(42)
n_simulations = 100_000

# Site visitors is symmetric (regular)
traffic_dist = stats.norm(loc=50000, scale=5000)

# Conversion charge is a likelihood bounded between 0 and 1 (beta)
# Imply = 5 / (5 + 95) = 5%
conversion_dist = stats.beta(a=5, b=95)

# Common order worth is optimistic and right-skewed (gamma)
# Imply = 15 * 4 = $60
aov_dist = stats.gamma(a=15, scale=4)

# 2. Draw N samples
traffic_samples = traffic_dist.rvs(measurement=n_simulations)
conversion_samples = conversion_dist.rvs(measurement=n_simulations)
aov_samples = aov_dist.rvs(measurement=n_simulations)

# 3. Vectorized propagation by means of the enterprise logic
revenue_samples = traffic_samples * conversion_samples * aov_samples

# 4. Extract threat percentiles
mean_rev = np.imply(revenue_samples)

# 5% likelihood of constructing lower than this
p5_rev = np.percentile(revenue_samples, 5)

# 5% likelihood of constructing greater than this
p95_rev = np.percentile(revenue_samples, 95)

print(f"Probabilistic Imply Income:  ${mean_rev:,.2f}")
print(f"fifth Percentile (Draw back):    ${p5_rev:,.2f}")
print(f"ninety fifth Percentile (Upside):     ${p95_rev:,.2f}")

Ausgabe:

Probabilistic Imply Income:  $150,153.16
fifth Percentile (Draw back):    $51,294.37
ninety fifth Percentile (Upside):     $300,734.30

Während der durchschnittliche Umsatz unserer ursprünglichen Punktschätzung entspricht (ca. 150.000 US-Greenback), zeigt die Monte-Carlo-Simulation, dass die Wahrscheinlichkeit, dass unser Umsatz darunter liegt, bei 5 % liegt 97.180 $ aufgrund der gemeinsamen Volatilität von Site visitors, Conversion und Bestellwert. Punktschätzungen verbergen diese betriebliche Gefährdung vollständig.

# 3. Sensitivitätsanalyse über Parameter-Sweeps

Bei der Szenarioanalyse möchten Sie möglicherweise herausfinden, wie empfindlich Ihr Abwärtsrisiko auf Änderungen bestimmter Inputannahmen reagiert. Wenn Sie beispielsweise die Kundenakquisekosten (CAC) modellieren, möchten Sie verstehen, wie eine Verschiebung unserer Marketingvolatilität (Standardabweichung) unseren Worst-Case-CAC (95. Perzentil) verschiebt.

Während Sie für jede Konfiguration eine vollständige Monte-Carlo-Simulation durchführen könnten, scipy.stats bietet eine viel schnellere, mathematisch elegante Abkürzung: die Prozentpunktfunktion (.ppf()), die die Umkehrung der kumulativen Verteilungsfunktion (CDF) ist.

Indem wir unsere Parameter durchsuchen, die Verteilungen einfrieren und ausführen .ppf()können wir die genauen Perzentilgrenzen analytisch in Mikrosekunden berechnen, ohne eine einzige Zufallsstichprobe zu generieren.

Die Simulation tausender Stichproben für jede Parameterpermutation, nur um ein Perzentil zu finden, ist höchst ineffizient und rechenintensiv.

import numpy as np
import scipy.stats as stats

mean_cac = 50.0
volatilities = (5.0, 10.0, 15.0, 20.0)

# Operating an enormous random simulation simply to extract a single percentile
for vol in volatilities:
    samples = stats.norm.rvs(loc=mean_cac, scale=vol, measurement=100000)
    p95_clunky = np.percentile(samples, 95)
    print(f"Std: {vol} -> ninety fifth Percentile CAC: ${p95_clunky:.2f} (sampled)")

Beispielausgabe:

Std: 5.0 -> ninety fifth Percentile CAC: $58.23 (sampled)
Std: 10.0 -> ninety fifth Percentile CAC: $66.53 (sampled)
Std: 15.0 -> ninety fifth Percentile CAC: $74.65 (sampled)
Std: 20.0 -> ninety fifth Percentile CAC: $82.85 (sampled)

Durch die Nutzung der .ppf() Mithilfe der Methode für unsere eingefrorenen Verteilungen berechnen wir sofort den genauen analytischen Schwellenwert.

import scipy.stats as stats

mean_cac = 50.0
volatilities = (5.0, 10.0, 15.0, 20.0)

# The SciPy Means: Sweep parameters and compute actual percentiles utilizing .ppf()
for vol in volatilities:
    # 1. Freeze the distribution
    cac_dist = stats.norm(loc=mean_cac, scale=vol)

    # 2. Compute actual ninety fifth percentile analytically
    p95_exact = cac_dist.ppf(0.95)

    # 3. Compute the likelihood of exceeding an excessive threshold of $80
    p_exceed_80 = cac_dist.sf(80.0)

    print(f"Volatility (Std): ${vol:02.0f} | ninety fifth Percentile CAC: ${p95_exact:.2f} | P(CAC > $80): {p_exceed_80:.2%}")

Ausgabe:

Volatility (Std): $05 | ninety fifth Percentile CAC: $58.22 | P(CAC > $80): 0.00%
Volatility (Std): $10 | ninety fifth Percentile CAC: $66.45 | P(CAC > $80): 0.13%
Volatility (Std): $15 | ninety fifth Percentile CAC: $74.67 | P(CAC > $80): 2.28%
Volatility (Std): $20 | ninety fifth Percentile CAC: $82.90 | P(CAC > $80): 6.68%

Dieser Durchlauf legt unsere Grenzen offen: Wenn unsere Akquisitionsvolatilität von 5 US-Greenback auf 20 US-Greenback steigt, steigen unsere Akquisitionskosten im 95. Perzentil 58,22 $ Zu 82,90 $und das Risiko, unser maximales Akquisitionsbudget von 80 US-Greenback zu überschreiten, steigt 0,00 % Zu 6,68 %.

# 4. Modellierung des Tail-Risikos mit Heavy-Tailed-Verteilungen

Ein häufiger Fehler bei der Szenarioplanung ist die Annahme, dass jeder kontinuierliche Datensatz einer Customary-Gauß-Verteilung (Normalverteilung) folgt. Obwohl die Normalverteilung einfach zu modellieren ist, weist sie extrem dünne Enden auf. In der realen Welt kommt es zu Systemausfällen, finanziellen Verlusten und API-Latenzen schwerschwänzig: Extremereignisse treten weitaus häufiger auf, als ein Gaußsches Modell jemals vorhersagen würde.

Um unsere Systeme ordnungsgemäß auf das Extremrisiko zu testen, müssen wir naive Normalannahmen durch stark ausgeprägte Alternativen wie Scholar’s t (stats.t), logarithmisch regular (stats.lognorm) oder Pareto (stats.pareto).

Mit der .match() Methode in scipy.statskönnen wir sowohl Regular- als auch Heavy-Tail-Verteilungen an historische Beobachtungen anpassen und dann die Überlebensfunktion verwenden (.sf()), um die realistische Wahrscheinlichkeit von Ausfällen im schlimmsten Fall zu quantifizieren.

Wenn Sie mit stark ausgeprägten Daten arbeiten, werden Sie durch die Modellierung mit einer Normalverteilung völlig blind gegenüber extremen Ereignissen nach unten:

import numpy as np
import scipy.stats as stats

# Artificial historic latency knowledge with heavy tails (Scholar's t with df=3)
np.random.seed(42)
latencies = stats.t(df=3, loc=200, scale=40).rvs(measurement=1000)

# Naively assuming regular distribution with out testing match
mean_lat, std_lat = np.imply(latencies), np.std(latencies)
prob_outage_clunky = 1 - stats.norm.cdf(400, loc=mean_lat, scale=std_lat)
print(f"Naive Gaussian P(Latency > 400ms): {prob_outage_clunky:.6%}")

Ausgabe:

Naive Gaussian P(Latency > 400ms): 0.046321%

Indem wir eine Scholar-t-Verteilung an die Normalverteilung anpassen, können wir die Güte der Anpassung explizit vergleichen und Extremrisiken genau berechnen.

import numpy as np
import scipy.stats as stats

# Artificial historic latency knowledge (heavy-tailed)
np.random.seed(42)
latencies = stats.t(df=3, loc=200, scale=40).rvs(measurement=1000)

# 1. Match regular and heavy-tailed Scholar's t distributions to historic knowledge
norm_params = stats.norm.match(latencies)
t_params = stats.t.match(latencies)

# 2. Freeze the fitted fashions
fitted_norm = stats.norm(*norm_params)
fitted_t = stats.t(*t_params)

# 3. Calculate likelihood of exceeding a extreme timeout threshold of 400ms
prob_outage_norm = fitted_norm.sf(400)
prob_outage_t = fitted_t.sf(400)

# 4. Calculate the 99th percentile response time for SLA stress-testing
p99_norm = fitted_norm.ppf(0.99)
p99_t = fitted_t.ppf(0.99)

print(f"Regular SLA (99th percentile):   {p99_norm:.2f} ms | P(Latency > 400ms): {prob_outage_norm:.4%}")
print(f"Heavy-t SLA (99th percentile):  {p99_t:.2f} ms | P(Latency > 400ms): {prob_outage_t:.4%}")

Ausgabe:

Regular SLA (99th percentile):   340.14 ms | P(Latency > 400ms): 0.0463%
Heavy-t SLA (99th percentile):  368.97 ms | P(Latency > 400ms): 0.6161%

Der Unterschied ist spürbar. Unter der naiven Gaußschen Annahme wird vorhergesagt, dass ein Ausfall mit hoher Latenz von mehr als 400 ms ein seltenes Ereignis ist (vorkommt). 0,0463 % der Zeit). In Wirklichkeit zeigt das Heavy-Tailed-Modell, dass die Ausfallwahrscheinlichkeit gleich ist 0,6161 % – über 13x häufiger. Durch die Anpassung von Distributionen mit starken Schwankungen werden Sie nicht von scheinbar „seltenen“ Fehlern überrascht.

# 5. Bootstrapping-Konfidenzintervalle für Szenariometriken

Wenn Sie eine Simulation ausführen oder einen kleinen historischen Datensatz analysieren, berechnen Sie eine zusammenfassende Metrik (z. B. das 90. Perzentil der Auftragsgrößen oder die mittleren Betriebskosten). Aber wie stabil ist diese Metrik? Was wäre, wenn Ihre historische Stichprobe etwas anders wäre?

Die analytische Berechnung von Konfidenzintervallen für nicht standardmäßige Metriken (wie ein Verhältnis oder ein Perzentil) ist mathematisch komplex. In der Vergangenheit haben Praktiker klobige verschachtelte Schleifen geschrieben, um Daten manuell erneut abzutasten.

SciPy 1.7 führte den neuesten Stand der Technik ein scipy.stats.bootstrap Funktion. In einem einzigen, hochoptimierten Funktionsaufruf können Sie robuste, nicht parametrische, voreingenommene und beschleunigte (BCa) Konfidenzintervalle für berechnen beliebig willkürliche zusammenfassende Statistik, ohne die Annahme einer zugrunde liegenden Verteilung.

Das Schreiben verschachtelter Schleifen zum erneuten Abtasten, Berechnen von Statistiken und Ermitteln der Quantile Ihrer Bootstrap-Verteilung ist langsam, ausführlich und kann Bias- oder Beschleunigungsanpassungen nicht verarbeiten.

import numpy as np

# A small pattern of fifty transaction quantities
np.random.seed(42)
transactions = np.random.lognormal(imply=4, sigma=0.8, measurement=50)

# Guide bootstrap loop
boot_medians = ()
for _ in vary(10000):
    pattern = np.random.alternative(transactions, measurement=len(transactions), change=True)
    boot_medians.append(np.median(pattern))

ci_low = np.percentile(boot_medians, 2.5)
ci_high = np.percentile(boot_medians, 97.5)

print(f"Guide Bootstrap Median CI: ({ci_low:.2f}, {ci_high:.2f})")

Ausgabe:

Guide Bootstrap Median CI: (36.47, 60.12)

Im Gegensatz dazu übergeben wir unsere Daten- und Statistikfunktion direkt an stats.bootstrapSciPy berechnet ein Bias-korrigiertes Konfidenzintervall in Millisekunden.

import numpy as np
import scipy.stats as stats

# A small pattern of fifty transaction quantities
np.random.seed(42)
transactions = np.random.lognormal(imply=4, sigma=0.8, measurement=50)

# Outline the statistic we wish to assemble a CI for (should settle for knowledge as a sequence)
def my_median(data_seq):
    return np.median(data_seq)

# Execute stats.bootstrap utilizing BCa technique, passing knowledge as a tuple containing our array
bootstrap_res = stats.bootstrap(
    (transactions,),
    statistic=my_median,
    n_resamples=9999,
    confidence_level=0.95,
    technique='BCa'
)

ci = bootstrap_res.confidence_interval

print(f"Pattern Median transaction measurement: ${np.median(transactions):.2f}")
print(f"95% Confidence Interval (BCa):  (${ci.low:.2f}, ${ci.excessive:.2f})")
print(f"Customary Error of the Median:   ${bootstrap_res.standard_error:.4f}")

Ausgabe:

95% Confidence Interval (BCa):  ($36.47, $60.12)
Customary Error of the Median:   $5.8819

Beachten Sie, dass die BCa-Methode im Vergleich zum naiven manuellen Bootstrap ein engeres und äußerst genaues, mathematisch fundiertes Konfidenzintervall zurückgab. BCa passt sich automatisch an Schiefe und Verzerrungen im zugrunde liegenden Datensatz an und liefert eine prinzipielle Unsicherheit, die über jedem Szenario oder jeder Stichprobenschätzung liegt.

# Zusammenfassung

Der Übergang vom einfachen Punktschätzungsdenken zum ausgereiften Verteilungsdenken ist einer der wirkungsvollsten Schritte, die Sie als Datenwissenschaftler unternehmen können.

Durch die Einbeziehung dieser fünf scipy.stats Wenn Sie Methods in Ihren Modellierungsworkflow integrieren, können Sie äußerst belastbare, mathematisch fundierte Szenariosysteme entwerfen:

• Das Einfrieren von Verteilungen kapselt Ihre Geschäftsannahmen in saubere, austauschbare Szenarioobjekte
• Monte-Carlo-Probenahme über .rvs() propagiert multivariable Unsicherheiten nahtlos in vektorisierten Hochgeschwindigkeits-C-Erweiterungen
• Parameter-Sweeps mit .ppf() Berechnen Sie präzise Risikoschwellenwerte und Perzentile analytisch in Mikrosekunden
• Schwerer Beschlag mit .match() Und .sf() schützt Ihre Produktionsarchitekturen vor katastrophalen Black-Swan-Ereignissen
• BCa Bootstrapping mit stats.bootstrap Platziert robuste, verteilungsfreie Konfidenzgrenzen über jeder Szenariometrik

Wenn Sie das nächste Mal damit beauftragt werden, Systeme einem Stresstest zu unterziehen oder Geschäftsrisiken unter Unsicherheit abzuschätzen, überspringen Sie die komplexen externen Abhängigkeiten. Schnappen Sie sich Ihren standardmäßigen wissenschaftlichen Python-Stack und nutzen Sie die Leistungsfähigkeit von scipy.statsund lassen Sie die Simulationen laufen!

Matthew Mayo (@mattmayo13) hat einen Grasp-Abschluss in Informatik und ein Diplom in Information Mining. Als geschäftsführender Herausgeber von KDnuggets & Statistikund Mitherausgeber bei Beherrschung des maschinellen LernensZiel von Matthew ist es, komplexe datenwissenschaftliche Konzepte zugänglich zu machen. Zu seinen beruflichen Interessen zählen die Verarbeitung natürlicher Sprache, Sprachmodelle, Algorithmen für maschinelles Lernen und die Erforschung neuer KI. Seine Mission ist es, das Wissen in der Datenwissenschaftsgemeinschaft zu demokratisieren. Matthew programmiert seit seinem sechsten Lebensjahr.