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

import numpy as np
import math
 
def g():
 
    g1=0
    for i in range(number):
        g1=g1+(y[i]-e[4][i])*(y[i]-e[4][i])
 
    return g1
 
def f(ax):
 
    for n in range(number_eq):
        for i in range(number):
            f1=ax[n][0] 
            for j in range(number_var):
                n1=link[n][j+1]
                f1=f1+ax[n][j+1]*e[n1][i]
            n1=link[n][0]
            e[n1][i]=1/(1+math.exp(-f1)) 
 
number=100
number_eq=3
number_var=2
 
x=[]
for i in range(100):
   x.append([2*i/200,(2*i+1)/200])
 
y=[]
for i in range(50):
   y.append(1)
for i in range(50):
   y.append(0)
 
 
link=[[2,0,1],[3,0,1],[4,2,3]]
e=np.zeros([5,number])
a=np.zeros([number_eq,number_var+1])
a[0][0]=0.4
a[0][1]=0.5
a[0][2]=0.7
a[1][0]=0.6
a[1][1]=0.8
a[1][2]=0.2
a[2][0]=0.4
a[2][1]=0.5
a[2][2]=0.7
 
 
for i in range(number):
    e[0][i]=x[i][0]
    e[1][i]=x[i][1]
 
ax=np.zeros([number_eq,number_var+1])
df=np.zeros([number_eq,number_var+1])
 
tr=1
while tr<500:
    tr=tr+1
    f(a)
    g1=g()
 
    for eq1 in range(number_eq): 
      for m in range(number_var+1):
 
 
 
        for eq2 in range(number_eq):
            for i in range(number_var+1):
                ax[eq2][i]=a[eq2][i]
        h=0.01
        ax[eq1][m]=ax[eq1][m]+h
        f(ax)
        g2=g()
        df[eq1][m]=(g2-g1)/h
 
    sig=0
    for eq1 in range(number_eq): 
      for m in range(number_var+1):
        sig=sig+df[eq1][m]*df[eq1][m]
 
    for eq1 in range(number_eq): 
      for m in range(number_var+1):
        ax[eq1][m]=a[eq1][m]-0.01*g1*df[eq1][m]/sig
 
    f(ax)
    g2=g()
 
    if g2<g1:
      for eq1 in range(number_eq): 
        for m in range(number_var+1):
            a[eq1][m]=ax[eq1][m]
    if g2>g1:tr=1000
 
for i in range(number):
    z1=a[0][0]+a[0][1]*x[i][0]+a[0][2]*x[i][1]
    e1=0
    if z1>0:e1=1
    z2=a[1][0]+a[1][1]*x[i][0]+a[1][2]*x[i][1]
    e2=0
    if z2>0:e2=1
    z3=a[2][0]+a[2][1]*e1+a[2][2]*e2
    e3=0
    if z3>0:e3=1
    print(str(e3)+","+str(y[i]))