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

#*****************************************************************************80
#
## ford_normal uses the normal equations to model Ford 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
  from ford_data import ford_data

  print ( '' )
  print ( 'ford_normal:' )
  print ( '  Given mileage x and price y for used Fords,' )
  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.
#
  x, y = ford_data ( )
  n = len ( x )
#
#  Create the matrix A:
#
  A = np.zeros ( [ n, 2 ] )
  A[:,0] = 1;
  A[:,1] = x;
  print ( A )
#
#  Create the system matrix A'A:
#
  AtA = np.matmul ( A.transpose(), A )
  print ( AtA )
#
#  Create the system right hand side A' y:
#
  Aty = np.matmul ( A.transpose(), y )
#
#  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 = normalized mileage' )
  print ( '  y = normalized price' )
  print ( '  Seek relationship y = b + m * x' )
  print ( '  b = %g' % ( b ) )
  print ( '  m = %g' % ( m ) )
  e = np.sum ( ( y - b - m * x )**2 )
  print ( '  Sum-of-squares error is %g' % ( e ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'ford_normal:' )
  print ( '  Normal end of execution.' )
  return

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