OP有两个问题。
第一个问题是代码试图用正态分布来拟合随机样本。这是错误的。但是，可以拟合直方图的输出。如下面的代码所示。最好使用scipy.stats.norm.fit()，它允许拟合随机样本。这也显示了。
第二个问题是sigma-shape。这里curve_fit实际上预期y数据上的误差，这自然需要y数据的形状。应该做的是：提供拟合的起始值。这也如下所示。
代码如下所示：

import matplotlib.pyplot as plt
import numpy as np

from scipy.stats import norm
from scipy.optimize import curve_fit

mean_o = list()
sigma_o = list()
y3 = list()

### generate some data

for i in range( 5 ):
    y3.append( norm.rvs( size=150000 ) )
y3 = np.transpose( y3 )

for i in range(5):
    mean_o.append( np.mean( y3[ :, i ] ) )
    sigma_o.append(  np.std( y3[ :, i ] ) )

## Histograms

# Number of bins

Nbins=100
binwidth = np.zeros(5)

# Fitting curves

def gaussian( x, a , mean, sigma ):
  return a * np.exp( -( ( x - mean )**2 / (2 * sigma**2 ) ) )

fig = plt.figure()
ax = { i : fig.add_subplot( 2, 3, i+1) for i in range( 5 ) }

for i in range(5):
    ymin = min(y3[:,i])
    ymax = max(y3[:,i])
    binwidth[i] = ( ymax - ymin) / Nbins
    bins_plot = np.arange( ymin, ymax + binwidth[i], binwidth[i])
    histdata = ax[i].hist(
        y3[:,i],
        bins=bins_plot,
        label='bin '+str(i)
    )

    range_fit = np.linspace( ymin, ymax, 250)
    # Fitting and plot version 1
    popt, pcov = curve_fit(
        gaussian,
        0.5 * ( histdata[1][:-1] + histdata[1][1:] ),
        histdata[0],
        p0=( max(histdata[0]), mean_o[i], sigma_o[i] ) )
    ax[i].plot(range_fit, gaussian( range_fit, *popt ) )
    ax[i].axvline( x=mean_o[i], ls=':', c='r' )

    # Fitting and plot version 2
    params = norm.fit( y3[ ::, i ], loc=mean_o[i], scale=sigma_o[i] )
    nth = gaussian(
        range_fit,
        len( y3[::, i]) * binwidth[i] / np.sqrt( 2 * np.pi ),
        *params
    )
    ax[i].plot(range_fit, nth, ls="--" )
plt.tight_layout()
plt.show()