Anpassung einer Weibull-Verteilung mit Scipy

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). 

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.  Prost!

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

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:

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()

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()

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:

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()

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.