#! /usr/bin/env python3
#
def china_normal ( ):

#*****************************************************************************80
#
## china_normal uses the normal equations to model some income data.
#
#  Discussion:
#
#    Normal equations to solve for optimal (b,m) to fit (x,y) data
#    y = b + m*x
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    09 October 2019
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  print ( '' )
  print ( 'china_normal:' )
  print ( '  Given year x and income y,' )
  print ( '  seek (m,b) so that y=b+mx approximates the data.' )
  print ( '  Use the normal equations to solve for best b and m.' )
#
#  Get the (normalized) Ford data.
#
  data = np.loadtxt ( 'china_data.txt' )
  year = data[:,0]
  income = data[:,1]
  n = len ( year )
#
#  Create the matrix A:
#
  A = np.zeros ( [ n, 2 ] )
  A[:,0] = 1;
  A[:,1] = year;
#
#  Create the system matrix A'A:
#
  AtA = np.matmul ( A.transpose(), A )
#
#  Create the system right hand side A' y:
#
  Aty = np.matmul ( A.transpose(), income )
#
#  Solve the linear system AtA * bm = Aty
#
  bm = np.linalg.solve ( AtA, Aty )
  b = bm[0]
  m = bm[1]
#
#  Report the results.
#
  print ( '' )
  print ( '  Estimating model parameters:' )
  print ( '  x = year' )
  print ( '  y = income' )
  print ( '  Seek relationship y = b + m * x' )
  print ( '  b = %g' % ( b ) )
  print ( '  m = %g' % ( m ) )
  e = np.sum ( ( income - b - m * year )**2 )
  print ( '  Sum-of-squares error is %g' % ( e ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'china_normal:' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  china_normal ( )