#NEWLINE# @ skip_if_32bit #NEWLINE# def test_incremental_variance_numerical_stability( ) : #NEWLINE# #INDENT# def np_var( A ) : #NEWLINE# #INDENT# return A . var ( axis = 0 ) #NEWLINE# #UNINDENT# def one_pass_var( X ) : #NEWLINE# #INDENT# n = X . shape [ 0 ] #NEWLINE# exp_x2 = ( ( X ** 2 ) . sum ( axis = 0 ) / n ) #NEWLINE# expx_2 = ( ( X . sum ( axis = 0 ) / n ) ** 2 ) #NEWLINE# return ( exp_x2 - expx_2 ) #NEWLINE# #UNINDENT# def two_pass_var( X ) : #NEWLINE# #INDENT# mean = X . mean ( axis = 0 ) #NEWLINE# Y = X . copy ( ) #NEWLINE# return np . mean ( ( ( Y - mean ) ** 2 ) , axis = 0 ) #NEWLINE# #UNINDENT# def naive_mean_variance_update( x , last_mean , last_variance , last_sample_count ) : #NEWLINE# #INDENT# updated_sample_count = ( last_sample_count + 1 ) #NEWLINE# samples_ratio = ( last_sample_count / float ( updated_sample_count ) ) #NEWLINE# updated_mean = ( ( x / updated_sample_count ) + ( last_mean * samples_ratio ) ) #NEWLINE# updated_variance = ( ( last_variance * samples_ratio ) + ( ( ( x - last_mean ) * ( x - updated_mean ) ) / updated_sample_count ) ) #NEWLINE# return ( updated_mean , updated_variance , updated_sample_count ) #NEWLINE# #UNINDENT# tol = 200 #NEWLINE# n_features = 2 #NEWLINE# n_samples = 10000 #NEWLINE# x1 = np . array ( 100000000.0 , dtype = np . float64 ) #NEWLINE# x2 = np . log ( 1e-05 , dtype = np . float64 ) #NEWLINE# A0 = ( x1 * np . ones ( ( ( n_samples // 2 ) , n_features ) , dtype = np . float64 ) ) #NEWLINE# A1 = ( x2 * np . ones ( ( ( n_samples // 2 ) , n_features ) , dtype = np . float64 ) ) #NEWLINE# A = np . vstack ( ( A0 , A1 ) ) #NEWLINE# if ( np . abs ( ( np_var ( A ) - two_pass_var ( A ) ) ) . max ( ) < 1e-06 ) : #NEWLINE# #INDENT# stable_var = np_var #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# stable_var = two_pass_var #NEWLINE# #UNINDENT# assert_greater ( np . abs ( ( stable_var ( A ) - one_pass_var ( A ) ) ) . max ( ) , tol ) #NEWLINE# ( mean , var , n ) = ( A0 [ 0 , : ] , np . zeros ( n_features ) , ( n_samples // 2 ) ) #NEWLINE# for i in range ( A1 . shape [ 0 ] ) : #NEWLINE# #INDENT# ( mean , var , n ) = naive_mean_variance_update ( A1 [ i , : ] , mean , var , n ) #NEWLINE# #UNINDENT# assert_equal ( n , A . shape [ 0 ] ) #NEWLINE# assert_greater ( np . abs ( ( A . mean ( axis = 0 ) - mean ) ) . max ( ) , 1e-06 ) #NEWLINE# assert_greater ( np . abs ( ( stable_var ( A ) - var ) ) . max ( ) , tol ) #NEWLINE# ( mean , var , n ) = ( A0 [ 0 , : ] , np . zeros ( n_features ) , ( n_samples // 2 ) ) #NEWLINE# for i in range ( A1 . shape [ 0 ] ) : #NEWLINE# #INDENT# ( mean , var , n ) = _incremental_mean_and_var ( A1 [ i , : ] . reshape ( ( 1 , A1 . shape [ 1 ] ) ) , mean , var , n ) #NEWLINE# #UNINDENT# assert_equal ( n , A . shape [ 0 ] ) #NEWLINE# assert_array_almost_equal ( A . mean ( axis = 0 ) , mean ) #NEWLINE# assert_greater ( tol , np . abs ( ( stable_var ( A ) - var ) ) . max ( ) )