ん６６６６ - 経済統計練習帳

import numpy as np
import math
 
def df(a1,m):
 
    f1=f(a1)
 
    h=0.01
    a1[m]=a1[m]+h
    f2=f(a1)
 
    return (f2-f1)/h
 
 
 
def f(a1):
 
    e=[]
    for tr in range(5): 
        z=a[0]
        for i in range(2):
            z=z+a1[i+1]*x[tr][i]
        e.append(z)
 
    sig=0
    for tr in range(5):
        sig=sig+(y[tr]-1/(1+math.exp(-e[tr])))*(y[tr]-1/(1+math.exp(-e[tr])))
 
    return sig
 
a=[0.4,0.5,0.2]
x=[[0.7,0.5],[0.4,0.6],[0.23,0.4],[0.4,0.6],[0.9,0.1]]
y=[0,0,0,1,1]
 
trs=0
while trs<100:
    trs=trs+1
    dg=[0,0,0]
 
    for m in range(3):
        dg[m]=df(a,m)
 
    e=-0.01*f(a)
 
    mu=0
    for m in range(3):
        mu=mu+dg[m]*dg[m]
    mu=math.sqrt(mu)
 
    delta=[0,0,0]
    for m in range(3):
        delta[m]=e*dg[m]/mu
 
    ax=[0,0,0]
    for m in range(3):
        ax[m]=a[m]+delta[m]
 
    f_de=f(a)
    f_delta=f(ax)
 
    if f_delta<f_de:
        for m in range(3):
           a[m]=ax[m]
 
    print(a)