「scipy.optimize.minimize」でパラメータ (atheta) を推定しました。私の手順は、最尤法を計算することと同じです。統計ソフトウェアが行うように、このパラメーターの標準誤差を計算したいと考えています。
パッケージ scikits.bootstrap を見つけましたが、カスタム関数の信頼区間を計算するのではなく、scipy 統計関数のみを計算するようです。
標準誤差を計算するにはどうすればよいですか?
これが私のコードです:
from __future__ import division
import numpy as np
import pandas
import scipy
from scipy import optimize
# import data
dir =
data =
#define function to minimize
def f(y, ns, vars):
atheta = y[:1]
tosum = 1/(np.exp(atheta)-np.exp(-atheta*vars))
sum = np.nansum(tosum,axis=1)
firstterm = tosum[:,[0]]
firsterm2 = firstterm.flatten()
lnp1 = np.log(firsterm2 * 1/sum)
return -np.sum(lnp1)
# this is the minimisation of the likelihood. It gives back atheta.
def main():
print '*'*80
print 'nouvelle execution'
print '*'*80
# data
ns = data['n'].values.astype('int')
vars = data.loc[:, ('R1', 'R2', 'R3', 'R4', 'R5', 'R6')].values
ns= np.array(ns, dtype=np.int)
vars= np.array(vars, dtype=np.float)
x0 = [-0.1]
result = scipy.optimize.minimize(f, x0, method = 'Nelder-Mead',
args = (ns, vars))
return result
if __name__ == "__main__":
print 'resultat du main = ', main()
データは次のようになります。
R1 R2 R3 R4 R5 R6 n
1 30.3 4.1 10.2 2.5 10.8 6
0.9 10.4 4.1 6.3 3.3 NaN 5
データには 25000 行があり、変数 R の数は最大 R24 になるため、これは単なるサンプルです。