こひｎ - 経済統計練習帳

import numpy as np
import math
 
def g():
 
    g3=0
    for i in range(number):
        g3=g3+(y[i]-e[0][i])*(y[i]-e[0][i])
 
    return g3
 
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):
   r1=np.random.random()
   r2=np.random.random()
   x.append([r1,r2])
 
y=[]
for i in range(100):
    r3=0
    h=0
    if 3*x[i][0]+2*x[i][1]>2.5:h=h+1  
    if h==1:r3=1
    y.append(r3)
 
link=[[1,3,4],[2,3,4],[0,1,2]]
 
e=np.zeros([5,number])
a=np.zeros([number_eq,number_var+1])
df=np.zeros([number_eq,number_var+1])
ax=np.zeros([number_eq,number_var+1])
 
for eq1 in range(number_eq): 
    for u in range(number_var+1):
        a[eq1][u]=np.random.random()
 
for i in range(number):
    e[3][i]=x[i][0]
    e[4][i]=x[i][1]
 
tr=0
while tr<100:
    tr=tr+1
 
    f(a)
    g1=g()
 
    for eqs in range(number_eq): 
        for mx in range(number_var+1):
 
 
 
            for eq1 in range(number_eq): 
                for m in range(number_var+1):
                    ax[eq1][m]=a[eq1][m]
 
            h=0.01
            ax[eqs][mx]=ax[eqs][mx]+h
            f(ax)
            g2=g()
            df[eqs][mx]=(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]
 
es=np.zeros([5,number])
for i in range(number):
    es[3][i]=x[i][0]
    es[4][i]=x[i][1]
 
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]*es[n1][i]
        n1=link[n][0]
        f2=0
        if f1>0:f2=1
        es[n1][i]=f2
 
c=0
for i in range(number):
    if y[i]==es[0][i]:c=c+1
 
print(c)