#! /usr/bin/env python # def truncated_normal_b_moments ( order_max, mu, sigma, b ): #*****************************************************************************80 # ## TRUNCATED_NORMAL_B_MOMENTS: moments of upper truncated Normal distribution. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 09 March 2015 # # Author: # # John Burkardt # # Reference: # # Phoebus Dhrymes, # Moments of Truncated Normal Distributions, # May 2005. # # Parameters: # # Input, integer ORDER_MAX, the maximum order of the moments. # 0 <= ORDER_MAX. # # Input, real MU, SIGMA, the mean and standard deviation of the # parent Normal distribution. # 0 < S. # # Input, real B, the upper truncation limit. # # Output, real VALUE[0:ORDER_MAX], the moment of the PDF. # import numpy as np value = np.zeros ( order_max + 1 ) for order in range ( 0, order_max + 1 ): value[order] = truncated_normal_b_moment ( order, mu, sigma, b ) return value def truncated_normal_b_moment ( order, mu, sigma, b ): #*****************************************************************************80 # ## TRUNCATED_NORMAL_B_MOMENT: a moment of upper truncated Normal distribution. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 09 March 2015 # # Author: # # John Burkardt # # Reference: # # Phoebus Dhrymes, # Moments of Truncated Normal Distributions, # May 2005. # # Parameters: # # Input, integer ORDER, the order of the moment. # 0 <= ORDER. # # Input, real MU, SIGMA, the mean and standard deviation of the # parent Normal distribution. # 0 < S. # # Input, real B, the upper truncation limit. # # Output, real VALUE, the moment of the PDF. # from normal_01_cdf import normal_01_cdf from normal_01_pdf import normal_01_pdf from r8_choose import r8_choose from sys import exit if ( order < 0 ): print ( '' ) print ( 'TRUNCATED_NORMAL_B_MOMENT - Fatal error!' ) print ( ' ORDER < 0.' ) exit ( 'TRUNCATED_NORMAL_B_MOMENT - Fatal error!' ) if ( sigma <= 0.0 ): print ( '' ) print ( 'TRUNCATED_NORMAL_B_MOMENT - Fatal error!' ) print ( ' SIGMA <= 0.0.' ) exit ( 'TRUNCATED_NORMAL_B_MOMENT - Fatal error!' ) b_h = ( b - mu ) / sigma b_pdf = normal_01_pdf ( b_h ) b_cdf = normal_01_cdf ( b_h ) if ( b_cdf == 0.0 ): print ( '' ) print ( 'TRUNCATED_NORMAL_B_MOMENT - Fatal error!' ) print ( ' PDF/CDF ratio fails, because B_CDF too small.' ) print ( ' B_PDF = %g' % ( b_pdf ) ) print ( ' B_CDF = %g' % ( b_cdf ) ) exit ( 'TRUNCATED_NORMAL_B_MOMENT - Fatal error!' ) f = b_pdf / b_cdf value = 0.0 irm2 = 0.0 irm1 = 0.0 for r in range ( 0, order + 1 ): if ( r == 0 ): ir = 1.0 elif ( r == 1 ): ir = - f else: ir = - b_h ** ( r - 1 ) * f + ( r - 1 ) * irm2 value = value + r8_choose ( order, r ) \ * mu ** ( order - r ) \ * sigma ** r * ir irm2 = irm1 irm1 = ir return value def truncated_normal_b_moment_test ( ): #*****************************************************************************80 # ## TRUNCATED_NORMAL_B_MOMENT_TEST tests TRUNCATED_NORMAL_B_MOMENT. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 09 March 2015 # # Author: # # John Burkardt # import numpy as np import platform test_num = 6 mu_test = np.array ( [ 0.0, 0.0, 0.0, 0.0, 0.0, 5.0 ] ) sigma_test = np.array ( [ 1.0, 1.0, 1.0, 2.0, 2.0, 1.0 ] ) b_test = np.array ( [ 0.0, 10.0, -10.0, 10.0, -10.0, 10.0 ] ) print ( '' ) print ( 'TRUNCATED_NORMAL_B_MOMENT_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' TRUNCATED_NORMAL_B_MOMENT evaluates moments' ) print ( ' of the upper Truncated Normal distribution.' ) for test in range ( 0, test_num ): mu = mu_test[test] sigma = sigma_test[test] b = b_test[test] print ( '' ) print ( ' Test = %d, Mu = %g, Sigma = %g, B = %g' \ % ( test, mu, sigma, b ) ) print ( ' Order Moment' ) print ( '\n' ) for order in range ( 0, 9 ): value = truncated_normal_b_moment ( order, mu, sigma, b ) print ( ' %2d %12g' % ( order, value ) ) # # Terminate. # print ( '' ) print ( 'TRUNCATED_NORMAL_B_MOMENT_TEST:' ) print ( ' Normal end of execution.' ) return if ( __name__ == '__main__' ): from timestamp import timestamp timestamp ( ) truncated_normal_b_moment_test ( ) timestamp ( )