Python (Scipy): Den Skalenparameter (Standardabweichung) einer Gaußschen Verteilung finden

Question

Es ist durchaus üblich, die Wahrscheinlichkeitsdichte eines Wertes innerhalb einer Wahrscheinlichkeitsdichtefunktion (PDF) zu berechnen. Stellen Sie sich vor wir haben eine Gaußsche Verteilung mit Mittelwert = 40, eine Standardabweichung von 5 und jetzt möchte die Wahrscheinlichkeitsdichte von Wert erhalten 32. Wir gehen würden wie:

In [1]: import scipy.stats as stats 
In [2]: print stats.norm.pdf(32, loc=40, scale=5) 
Out [2]: 0.022 

-> Die Wahrscheinlichkeitsdichte ist 2,2 %.

Aber jetzt, betrachten wir das umgekehrte Problem. Ich habe den Mittelwert, ich habe den Wert bei Probability-Dichte von 0,05 und ich möchte die Standardabweichung erhalten (d. H. Den Skalierungsparameter).

Was ich implementieren könnte, ist ein numerischer Ansatz: Erstellen Sie stats.norm.pdf mehrere Male mit dem Maßstab-Parameter schrittweise erhöht und nehmen Sie diese mit dem Ergebnis so nah wie möglich.

In meinem Fall gebe ich den Wert 30 als 5% -Markierung an. Also muss ich diese „Gleichung“ lösen:

stats.norm.pdf(30, loc=40, scale=X) = 0.05 

Es gibt eine scipy Funktion „PPF“ genannt, die die Umkehrung des PDF ist, so wird es den Wert für eine bestimmte Wahrscheinlichkeitsdichte zurück, aber ich haven‘ t hat eine Funktion gefunden, die den Skalierungsparameter zurückgibt.

Die Implementierung einer Iteration würde zu viel Zeit in Anspruch nehmen (sowohl beim Erstellen als auch beim Berechnen). Mein Skript wird riesig sein, also sollte ich Rechenzeit sparen. Könnte die Lambda-Funktion in diesem Fall helfen? Ich weiß ungefähr, was es tut, aber ich habe es bisher nicht benutzt. Irgendwelche Ideen dazu?

Vielen Dank!

Answer 1

Das wird zwei Lösungen sein, weil normale PDF symmetrisch um den Mittelwert ist. So wie es aussieht, müssen Sie eine Gleichung mit einer einzelnen Variablen lösen. Es wird keine geschlossene Lösung haben, so dass Sie z. scipy.optimize.fsolve um es zu lösen.

EDIT: siehe @ unutbu die Antwort für die geschlossene Form Lösung in Bezug auf Lambert W-Funktion.

Answer 2

Die normal probability density Funktion wird f von

gegeben

Bei f und x wir für lösen möchten. Lassen Sie uns sympy fragen, ob es die Gleichung lösen:

import sympy as sy 
from sympy.abc import x, y, sigma 

expr = (1/(sy.sqrt(2*sy.pi)*sigma) * sy.exp(-x**2/(2*sigma**2))) - y 
ans = sy.solve(expr, sigma)[0] 
print(ans) 
# sqrt(2)*exp(LambertW(-2*pi*x**2*y**2)/2)/(2*sqrt(pi)*y) 

So scheint es dort eine geschlossene gebildete Lösung in Bezug auf die LambertW function, W, die

z = W(z) * exp(W(z)) 

für alle komplexwertigen erfüllt z .

Wir sympy nutzen könnten, um auch für gegebene x und y das numerische Ergebnis zu finden, aber vielleicht wäre es schneller sein, um die numerische Arbeit zu tun mit scipy.special.lambertw:

import numpy as np 
import scipy.special as special 

def sigma_func(x, y): 
    results = set([np.real_if_close(
     np.sqrt(2)*np.exp(special.lambertw(-2*np.pi*x**2*y**2, k=k)/2) 
     /(2*np.sqrt(np.pi)*y)).item() for k in (0, -1)]) 
    results = [s for s in results if np.isreal(s)] 
    return results 

Im Allgemeinen ist die LambertW Funktion zurückkehrt komplexe Werte, aber wir sind nur interessiert an reellwertigen Lösungen für sigma. Per the docs, special.lambertw hat zwei teilweise echte Zweige, wenn k=0 und k=1. Der obige Code prüft, ob der zurückgegebene Wert (für diese zwei Zweige) real ist, und gibt eine Liste aller realen Lösungen zurück, falls sie existieren. Wenn keine echte Lösung existiert, wird dann eine leere Liste zurückgegeben. Das passiert, wenn der PDF-Wert y nicht für einen realen Wert von Sigma (für den angegebenen Wert von x) erreicht wird.

---

Sie können es wie folgt verwenden:

x = 30.0 
loc = 40.0 
y = 0.02 
s = sigma_func(loc-x, y) 
print(s) 
# [16.65817044316178, 6.830458938511113] 

import scipy.stats as stats 
for si in s: 
    assert np.allclose(stats.norm.pdf(x, loc=loc, scale=si), y) 

---

Im Beispiel Sie gab mit y = 0.025, gibt es keine Lösung für Sigma:

import numpy as np 
import scipy.stats as stats 
import matplotlib.pyplot as plt 

x = 30.0 
loc = 40.0 
y = 0.025 
s = np.linspace(5, 20, 100) 
plt.plot(s, stats.norm.pdf(x, loc=loc, scale=s)) 
plt.hlines(y, 4, 20, color='red') # the horizontal line y = 0.025 
plt.ylabel('pdf') 
plt.xlabel('sigma') 
plt.show() 

und so sigma_func(40-30, 0.025) gibt eine leere Liste:

In [93]: sigma_func(40-30, 0.025) 
Out [93]: [] 

---

Das Grundstück oben ist typisch in dem Sinne, dass, wenn y zu groß ist, gibt es Null Lösungen, bei dem Maximum der Kurve (nennen wir es y_max es) ist eine Lösung

In [199]: y_max = np.nextafter(np.sqrt(1/(np.exp(1)*2*np.pi*(10)**2)), -np.inf) 

In [200]: y_max 
Out[200]: 0.024197072451914336 

In [201]: sigma_func(40-30, y_max) 
Out[201]: [9.9999999776424] 

und y kleiner als die y_max gibt es zwei Lösungen.