webentwicklung-frage-antwort-db.com.de

Anpassung einer Weibull-Verteilung mit Scipy

Ich versuche, die Maximum-Likelihood-Verteilungsanpassung wiederherzustellen. Dies kann ich bereits in Matlab und R tun, aber jetzt möchte ich scipy verwenden. Insbesondere möchte ich die Weibull-Verteilungsparameter für meinen Datensatz schätzen.

Ich habe das versucht:

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

def weib(x,n,a):
    return (a / n) * (x / n)**(a - 1) * np.exp(-(x / n)**a)

data = np.loadtxt("stack_data.csv")

(loc, scale) = s.exponweib.fit_loc_scale(data, 1, 1)
print loc, scale

x = np.linspace(data.min(), data.max(), 1000)
plt.plot(x, weib(x, loc, scale))
plt.hist(data, data.max(), normed=True)
plt.show()

Und hol dir das:

(2.5827280639441961, 3.4955032285727947)

Und eine Distribution, die so aussieht:

Weibull distribution using Scipy

Ich habe das exponweib verwendet, nachdem ich dieses gelesen habe http://www.johndcook.com/distributions_scipy.html . Ich habe auch die anderen Weibull-Funktionen in scipy ausprobiert (nur für den Fall!).

In Matlab (mit dem Distribution Fitting Tool - siehe Screenshot) und in R (sowohl mit der MASS-Bibliotheksfunktion fitdistr als auch mit dem GAMLSS-Paket) erhalte ich a (loc) und b (scale) -Parameter wie 1.58463497 5.93030013. Ich glaube, dass alle drei Methoden die Maximum-Likelihood-Methode für die Verteilungsanpassung verwenden.

Weibull distribution using Matlab

Ich habe meine Daten gepostet hier wenn Sie es ausprobieren möchten! Und der Vollständigkeit halber verwende ich Python 2.7.5, Scipy 0.12.0, R 2.15.2 und Matlab 2012b.

Warum erhalte ich ein anderes Ergebnis?

32
kungphil

Ich vermute, dass Sie den Formparameter und die Skala der Weibull-Verteilung schätzen möchten, während Sie die Position festhalten. Bei der Korrektur von loc wird davon ausgegangen, dass die Werte Ihrer Daten und der Verteilung positiv sind und die Untergrenze bei Null liegt.

floc=0 hält den Standort auf Null, f0=1 behält den ersten Formparameter des exponentiellen Weibull bei eins.

>>> stats.exponweib.fit(data, floc=0, f0=1)
[1, 1.8553346917584836, 0, 6.8820748596850905]
>>> stats.weibull_min.fit(data, floc=0)
[1.8553346917584836, 0, 6.8820748596850549]

Die Passform im Vergleich zum Histogramm sieht ok aus, ist aber nicht sehr gut. Die Parameterschätzungen sind etwas höher als die von Ihnen erwähnten und stammen von R und matlab.

pdate

Das Diagramm, das jetzt zur Verfügung steht, kann ich am ehesten mit uneingeschränkter Anpassung, jedoch unter Verwendung von Startwerten, erreichen. Das Grundstück ist noch weniger hoch. Beachten Sie, dass Fit-Werte ohne vorangestelltes f als Startwerte verwendet werden.

>>> from scipy import stats
>>> import matplotlib.pyplot as plt
>>> plt.plot(data, stats.exponweib.pdf(data, *stats.exponweib.fit(data, 1, 1, scale=02, loc=0)))
>>> _ = plt.hist(data, bins=np.linspace(0, 16, 33), normed=True, alpha=0.5);
>>> plt.show()

exponweib fit

21
Josef

Es ist einfach zu überprüfen, welches Ergebnis das wahre MLE ist. Sie benötigen lediglich eine einfache Funktion, um die Wahrscheinlichkeit des Logs zu berechnen:

>>> def wb2LL(p, x): #log-likelihood
    return sum(log(stats.weibull_min.pdf(x, p[1], 0., p[0])))
>>> adata=loadtxt('/home/user/stack_data.csv')
>>> wb2LL(array([6.8820748596850905, 1.8553346917584836]), adata)
-8290.1227946678173
>>> wb2LL(array([5.93030013, 1.57463497]), adata)
-8410.3327470347667

Das Ergebnis von fit Methode von exponweib und R fitdistr (@Warren) ist besser und hat eine höhere Log-Wahrscheinlichkeit. Es ist wahrscheinlicher, dass es das wahre MLE ist. Es ist nicht überraschend, dass das Ergebnis von GAMLSS anders ist. Es ist ein völlig anderes statistisches Modell: das generalisierte additive Modell.

Immer noch nicht überzeugt? Wir können eine 2D-Konfidenzgrenzkurve um MLE zeichnen (siehe das Buch von Meeker und Escobar für Details). Multi-dimensional Confidence Region

Dies bestätigt erneut, dass array([6.8820748596850905, 1.8553346917584836]) die richtige Antwort ist, da die Wahrscheinlichkeit geringer ist als bei jedem anderen Punkt im Parameterraum. Hinweis:

>>> log(array([6.8820748596850905, 1.8553346917584836]))
array([ 1.92892018,  0.61806511])

BTW1, MLE-Anpassung scheint möglicherweise nicht genau mit dem Verteilungshistogramm übereinzustimmen. Eine einfache Möglichkeit, über MLE nachzudenken, ist, dass MLE die Parameterschätzung ist, die angesichts der beobachteten Daten am wahrscheinlichsten ist. Es muss nicht gut auf das Histogramm passen, da dies den mittleren quadratischen Fehler minimiert.

Übrigens scheinen Ihre Daten leptokurtisch und linksgerichtet zu sein, was bedeutet, dass die Weibull-Verteilung möglicherweise nicht gut zu Ihren Daten passt. Versuchen Sie es z. Gompertz-Logistic, die die Log-Wahrscheinlichkeit um weitere 100 erhöht. enter image description hereenter image description here Prost!

21
CT Zhu

Ich weiß, es ist ein alter Beitrag, aber ich hatte gerade ein ähnliches Problem und dieser Thread hat mir geholfen, es zu lösen. Dachte, meine Lösung könnte für andere wie mich hilfreich sein:

# Fit Weibull function, some explanation below
params = stats.exponweib.fit(data, floc=0, f0=1)
shape = params[1]
scale = params[3]
print 'shape:',shape
print 'scale:',scale

#### Plotting
# Histogram first
values,bins,hist = plt.hist(data,bins=51,range=(0,25),normed=True)
center = (bins[:-1] + bins[1:]) / 2.

# Using all params and the stats function
plt.plot(center,stats.exponweib.pdf(center,*params),lw=4,label='scipy')

# Using my own Weibull function as a check
def weibull(u,shape,scale):
    '''Weibull distribution for wind speed u with shape parameter k and scale parameter A'''
    return (shape / scale) * (u / scale)**(shape-1) * np.exp(-(u/scale)**shape)

plt.plot(center,weibull(center,shape,scale),label='Wind analysis',lw=2)
plt.legend()

Einige zusätzliche Informationen, die mir beim Verständnis geholfen haben:

Die Scipy Weibull-Funktion kann vier Eingabeparameter annehmen: (a, c), loc und scale. Wenn Sie die Position und den ersten Formparameter (a) festlegen möchten, wird dies mit floc = 0, f0 = 1 durchgeführt. Durch Anpassen erhalten Sie dann die Parameter c und scale, wobei c dem Formparameter der Zwei-Parameter-Weibull-Verteilung (häufig in der Winddatenanalyse verwendet) und scale dem Skalierungsfaktor entspricht.

Aus Dokumenten:

exponweib.pdf(x, a, c) =
    a * c * (1-exp(-x**c))**(a-1) * exp(-x**c)*x**(c-1)

Wenn a 1 ist, dann

exponweib.pdf(x, a, c) =
    c * (1-exp(-x**c))**(0) * exp(-x**c)*x**(c-1)
  = c * (1) * exp(-x**c)*x**(c-1)
  = c * x **(c-1) * exp(-x**c)

Daraus sollte der Zusammenhang mit der Weibull-Funktion „Windanalyse“ deutlicher hervorgehen

6
Peter9192

Ich war neugierig auf Ihre Frage und obwohl dies keine Antwort ist, vergleicht sie das Matlab -Ergebnis mit Ihrem Ergebnis und dem Ergebnis unter Verwendung von leastsq, das die beste Korrelation mit den angegebenen Daten ergab:

enter image description here

Der Code lautet wie folgt:

import scipy.stats as s
import numpy as np
import matplotlib.pyplot as plt
import numpy.random as mtrand
from scipy.integrate import quad
from scipy.optimize import leastsq

## my distribution (Inverse Normal with shape parameter mu=1.0)
def weib(x,n,a):
    return (a / n) * (x / n)**(a-1) * np.exp(-(x/n)**a)

def residuals(p,x,y):
    integral = quad( weib, 0, 16, args=(p[0],p[1]) )[0]
    penalization = abs(1.-integral)*100000
    return y - weib(x, p[0],p[1]) + penalization

#
data = np.loadtxt("stack_data.csv")


x = np.linspace(data.min(), data.max(), 100)
n, bins, patches = plt.hist(data,bins=x, normed=True)
binsm = (bins[1:]+bins[:-1])/2

popt, pcov = leastsq(func=residuals, x0=(1.,1.), args=(binsm,n))

loc, scale = 1.58463497, 5.93030013
plt.plot(binsm,n)
plt.plot(x, weib(x, loc, scale),
         label='weib matlab, loc=%1.3f, scale=%1.3f' % (loc, scale), lw=4.)
loc, scale = s.exponweib.fit_loc_scale(data, 1, 1)
plt.plot(x, weib(x, loc, scale),
         label='weib stack, loc=%1.3f, scale=%1.3f' % (loc, scale), lw=4.)
plt.plot(x, weib(x,*popt),
         label='weib leastsq, loc=%1.3f, scale=%1.3f' % Tuple(popt), lw=4.)

plt.legend(loc='upper right')
plt.show()
6

Hierauf gab es bereits hier und an anderen Stellen einige Antworten. likt in Weibull-Verteilung und die Daten in derselben Figur (mit numpy und scipy)

Es hat noch eine Weile gedauert, bis ich ein sauberes Spielzeugbeispiel gefunden hatte, daher wäre es nützlich, es zu posten.

from scipy import stats
import matplotlib.pyplot as plt

#input for pseudo data
N = 10000
Kappa_in = 1.8
Lambda_in = 10
a_in = 1
loc_in = 0 

#Generate data from given input
data = stats.exponweib.rvs(a=a_in,c=Kappa_in, loc=loc_in, scale=Lambda_in, size = N)

#The a and loc are fixed in the fit since it is standard to assume they are known
a_out, Kappa_out, loc_out, Lambda_out = stats.exponweib.fit(data, f0=a_in,floc=loc_in)

#Plot
bins = range(51)
fig = plt.figure() 
ax = fig.add_subplot(1, 1, 1)
ax.plot(bins, stats.exponweib.pdf(bins, a=a_out,c=Kappa_out,loc=loc_out,scale = Lambda_out))
ax.hist(data, bins = bins , normed=True, alpha=0.5)
ax.annotate("Shape: $k = %.2f$ \n Scale: $\lambda = %.2f$"%(Kappa_out,Lambda_out), xy=(0.7, 0.85), xycoords=ax.transAxes)
plt.show()
1
Keith

Ich hatte das gleiche Problem, fand aber diese Einstellung loc=0 im exponweib.fit hat die Pumpe für die Optimierung vorbereitet. Das war alles, was von @ user333700's answer benötigt wurde. Ich konnte Ihre Daten nicht laden - Ihr Datenlink verweist auf ein Bild, nicht auf Daten. Also habe ich stattdessen einen Test für meine Daten durchgeführt:

Plot of distribution fit to problematic (bimodal?) data

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

N=30
counts, bins = np.histogram(x, bins=N)
bin_width = bins[1]-bins[0]
total_count = float(sum(counts))

f, ax = plt.subplots(1, 1)
f.suptitle(query_uri)

ax.bar(bins[:-1]+bin_width/2., counts, align='center', width=.85*bin_width)
ax.grid('on')
def fit_pdf(x, name='lognorm', color='r'):
    dist = getattr(ss, name)  # params = shape, loc, scale
    # dist = ss.gamma  # 3 params

    params = dist.fit(x, loc=0)  # 1-day lag minimum for shipping
    y = dist.pdf(bins, *params)*total_count*bin_width
    sqerror_sum = np.log(sum(ci*(yi - ci)**2. for (ci, yi) in Zip(counts, y)))
    ax.plot(bins, y, color, lw=3, alpha=0.6, label='%s   err=%3.2f' % (name, sqerror_sum))
    return y

colors = ['r-', 'g-', 'r:', 'g:']

for name, color in Zip(['exponweib', 't', 'gamma'], colors): # 'lognorm', 'erlang', 'chi2', 'weibull_min', 
    y = fit_pdf(x, name=name, color=color)

ax.legend(loc='best', frameon=False)
plt.show()
1
hobs

die Reihenfolge von Ort und Maßstab ist im Code durcheinander gebracht:

plt.plot(x, weib(x, scale, loc))

der Parameter scale sollte an erster Stelle stehen.

0
Kaihua Cai