Adattamento di una distribuzione di Weibull utilizzando Scipy

Question

Sto tentando di ricreare la massima distribuzione di verosimiglianza, posso già farlo in Matlab e R, ma ora voglio usare scipy. In particolare, vorrei stimare i parametri di distribuzione di Weibull per il mio set di dati.

Ho provato questo:

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

E ottenere questo:

(2.5827280639441961, 3.4955032285727947) 

E una distribuzione che assomiglia a questo:

Sono stato con la exponweib dopo aver letto questo http://www.johndcook.com/distributions_scipy.html. Ho anche provato le altre funzioni di Weibull in Scipy (solo nel caso!).

In Matlab (utilizzando lo Strumento di distribuzione - vedi screenshot) e in R (utilizzando sia la funzione della libreria MASS fitdistr e il pacchetto GAMLSS) Ottengo i parametri (loc) eb (scala) più come 1.58463497 5.93030013. Credo che tutti e tre i metodi utilizzino il metodo della massima verosimiglianza per il fitting di distribuzione.

ho pubblicato i miei dati here se si desidera avere un andare! E per completezza sto usando Python 2.7.5, Scipy 0.12.0, R 2.15.2 e Matlab 2012b.

Perché sto ottenendo un risultato diverso !?

Answer 1

la mia ipotesi è che si vuole stimare il parametro di forma e la scala di distribuzione di Weibull, mantenendo la correzione posizione ed. La correzione loc presuppone che i valori dei dati e della distribuzione siano positivi con limite inferiore a zero.

floc=0 mantiene la posizione fissa su zero, f0=1 mantiene fisso il primo parametro di forma del weibull esponenziale.

>>> 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] 

La vestibilità rispetto all'istogramma sembra ok, ma non molto buona. Le stime dei parametri sono un po 'più alte di quelle che menzioni da R e MATLAB.

Aggiornamento

Il più vicino che posso ottenere per la trama che è ora disponibile sia con la misura illimitata, ma utilizzando i valori di partenza. La trama è ancora meno alta. I valori di nota in forma che non hanno un f in avanti vengono utilizzati come valori iniziali.

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

Answer 2

ero curioso di sapere la vostra domanda e, nonostante questo non è una risposta, si confronta il risultato Matlab con il risultato e con il risultato usando leastsq, che ha mostrato la migliore correlazione con i dati forniti:

il codice è il seguente:

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

Answer 3

È facile verificare che risultato è la vera MLE, basta una semplice funzione per calcolare log verosimiglianza:

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

Il risultato fit metodo di exponweib e R fitdistr (@Warren) è migliore e ha una maggiore probabilità di registro. È più probabile che sia il vero MLE. Non sorprende che il risultato di GAMLSS sia diverso. È un modello statistico completo e completo: Modello additivo generalizzato.

Ancora non convinto? Possiamo tracciare un diagramma del limite di confidenza 2D attorno a MLE, vedere il libro di Meeker e Escobar per i dettagli).

Anche in questo caso si verifica che array([6.8820748596850905, 1.8553346917584836]) è la risposta corretta poiché loglikelihood è inferiore a qualsiasi altro punto nello spazio dei parametri. Nota:

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

BTW1, MLE potrebbe non apparire correttamente per adattarsi all'istogramma di distribuzione. Un modo semplice per pensare a MLE è che MLE è la stima dei parametri più probabile dati i dati osservati. Non è necessario montare visivamente il pozzo dell'istogramma, questo sarà qualcosa che minimizza l'errore quadratico medio.

BTW2, i dati sembrano essere leptokurtic e left-skewed, il che significa che la distribuzione di Weibull potrebbe non adattarsi bene ai dati. Prova, ad es. Gompertz-Logistic, che migliora la probabilità di log da un altro circa 100. Cheers!

Answer 4

l'ordine di loc e la scala è incasinato nel codice:

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

il parametro di scala dovrebbe venire prima.

Answer 5

Ho avuto lo stesso problema, ma ho riscontrato che l'impostazione loc=0 in exponweib.fit inizializzava la pompa per l'ottimizzazione. Era tutto ciò che serviva da @ user333700's answer. Non è stato possibile caricare i dati: i tuoi punti data link su un'immagine, non sui dati. Così ho eseguito un test sui miei dati invece:

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

Answer 6

Lo so che è un vecchio post, ma ho appena affrontato un problema simile e questa discussione mi ha aiutato a risolverlo. Pensavo che la mia soluzione potrebbe essere utile per altri come me:

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

qualche informazione in più che mi ha aiutato a capire:

funzione SciPy Weibull può prendere quattro parametri di input: (a, c), loc e la scala. Si desidera correggere la loc e il primo parametro di forma (a), questo viene fatto con floc = 0, f0 = 1. L'adattamento fornirà quindi parametri c e scala, dove c corrisponde al parametro di forma della distribuzione a due parametri di Weibull (spesso utilizzata nell'analisi dei dati del vento) e la scala corrisponde al suo fattore di scala.

Da documenti:

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

Se a è 1, allora

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) 

Da questo, la relazione 'analisi vento' funzione di Weibull dovrebbe essere più chiaro

Answer 7

Ci devono sono state alcune risposte a questo già qui e in altri posti. likt in Weibull distribution and the data in the same figure (with numpy and scipy)

Mi ci è voluto ancora un po 'per trovare un esempio di giocattolo pulito, quindi penso che sarebbe utile postare.

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

Answer 8

Nella funzione forma, ci sono 3 parametri da considerare:

1. I parametri di forma: in questo caso, abbiamo due parametri di forma, che può essere fissato secondo f0 e f1 . (Provalo per te!). Di solito, il nome del parametro è indicato da f% d, dove d è il numero della forma.

2. Il parametro di posizione: utilizzare floc per risolvere questo problema. Se il valore di non corrisponde a di correzione, la media dei dati viene emessa come posizione.

3. Il parametro di scala: utilizzare fscale per risolvere questo problema.

Il ritorno di qualsiasi adattamento viene fuori in questo ordine.

seguito lungo le linee di @ Peter9192, ho trovato la soluzione migliore per un Weibull CDF di ~ 20-30 campioni di dati utilizzando la seguente: _,gamma,_alpha=scipy.stats.exponweib.fit(data,floc=0,f0=1)

La formula per il CDF è:

1-np.exp(-np.power(x/alpha,gamma)) I valori per dati Ho stimato utilizzando un metodo di stima KM, corrispondente a una distribuzione di Weibull mi ha dato buoni valori.

Per risolvere un come 1, I non ha trovato loc = 0, scale = 1 come il metodo migliore come si può vedere chiaramente nei 4 valori di parametro restituiti. In secondo luogo, l'utilizzo di gamma, alfa da esso non ha dato la media corretta di Weibull.

Infine, mi ha confermato che il metodo funziona meglio calcolando la media della distribuzione di Weibull utilizzando: Mean=alpha*scipy.special.gamma(1+(1/gamma)) I valori mi sono corrispondeva alla mia domanda.

È possibile controllare significa & formule CDF qui per riferimento: https://en.m.wikipedia.org/wiki/Weibull_distribution