Comment trouver la bonne courbe gaussienne avec des données?

Question

J'ai le code qui tire d'une gaussienne 1D:

import numpy as np 
from scipy.stats import norm 
from scipy.optimize import curve_fit 
import matplotlib.mlab as mlab 
import matplotlib.pyplot as plt 
import gauss 

# Beginning in one dimension: 
mean = 0; Var = 1; N = 1000 
scatter = np.random.normal(mean,np.sqrt(Var),N) 
scatter = np.sort(scatter) 
mu,sigma = norm.fit(scatter) 

Je mu et sigma obtenir l'aide norm.fit()

Maintenant, je voudrais obtenir mes paramètres en utilisant

xdata = np.linspace(-5,5,N) 
pop, pcov = curve_fit(gauss.gauss_1d,xdata,scatter) 

Le problème est que je ne sais pas comment mapper mes points épars (tirés d'un gaussien 1D) à la ligne x afin d'utiliser curve_fit. En outre, supposons que j'utilise simplement et mu et sigma comme précédemment.

Je intrigue en utilisant:

n, bins, patches = plt.hist(scatter,50,facecolor='green') 
y = 2*max(n)*mlab.normpdf(bins,mu,sigma) 
l = plt.plot(bins,y,'r--') 

plt.xlabel('x-coord') 
plt.ylabel('Occurrences') 
plt.grid(True) 
plt.show() 

Mais je dois deviner l'amplitude 2 * max (n). Cela fonctionne mais ce n'est pas robuste. Comment puis-je trouver l'amplitude sans deviner?

Answer 1

Pour éviter de deviner l'amplitude, appelez hist() avec normed=True, puis l'amplitude correspond à normpdf(). Pour faire un ajustement de courbe, je suggère d'utiliser non pas la densité mais la distribution cumulative: Chaque échantillon a une hauteur de 1/N, qui résume successivement jusqu'à 1. Ceci a l'avantage que vous n'avez pas besoin de grouper des échantillons dans des bacs.

import numpy as np 
from scipy.stats import norm 
from scipy.optimize import curve_fit 
import matplotlib.pyplot as plt 

# Beginning in one dimension: 
mean = 0; Var = 1; N = 100 
scatter = np.random.normal(mean,np.sqrt(Var),N) 
scatter = np.sort(scatter) 
mu1,sigma1 = norm.fit(scatter) # classical fit 

scat_sum = np.cumsum(np.ones(scatter.shape))/N # cumulative samples 
[mu2,sigma2],Cx = curve_fit(norm.cdf, scatter, scat_sum, p0=[0,1]) # curve fit 
print(u"norm.fit(): µ1= {:+.4f}, σ1={:.4f}".format(mu1, sigma1)) 
print(u"curve_fit(): µ2= {:+.4f}, σ2={:.4f}".format(mu2, sigma2)) 

fg = plt.figure(1); fg.clf() 
ax = fg.add_subplot(1, 1, 1) 
t = np.linspace(-4,4, 1000) 
ax.plot(t, norm.cdf(t, mu1, sigma1), alpha=.5, label="norm.fit()") 
ax.plot(t, norm.cdf(t, mu2, sigma2), alpha=.5, label="curve_fit()") 
ax.step(scatter, scat_sum, 'x-', where='post', alpha=.5, label="Samples") 
ax.legend(loc="best") 
ax.grid(True) 
ax.set_xlabel("$x$") 
ax.set_ylabel("Cumulative Probability Density") 
ax.set_title("Fit to Normal Distribution") 

fg.canvas.draw() 
plt.show() 

impressions

norm.fit(): µ1= +0.1534, σ1=1.0203 
curve_fit(): µ2= +0.1135, σ2=1.0444 

et parcelles