subroutine print_proje_blocks(hf_data, num_data, blocklength, filename)
implicit none
integer :: iunit, length, blocking_events, i
real(dp), intent(in) :: hf_data(:), num_data(:)
character(len=*), intent(in) :: filename
integer, intent(in) :: blocklength ! this is the minimum step size (e.g. 2)
real(dp) :: mean1, error1, eie1, mean2, error2, eie2
real(dp) :: covariance, mean_proje, final_error, final_eie
real(dp), allocatable :: this(:), that(:) ! stores this array
iunit = get_free_unit()
open(iunit, file=filename)
write(iunit, "(A)") "# Block_length Mean Error Error_in_Error"
length = size(hf_data, 1)
allocate(this(length))
allocate(that(length))
that = hf_data !Denominator data
this = num_data !Numerator data
call analyze_data(this, mean1, error1, eie1) !means & error in num
call analyze_data(that, mean2, error2, eie2) !means & error in denom
covariance = calc_covariance(that, this)
mean_proje = mean1 / mean2
final_error = abs(mean_proje) * &
sqrt((error2 / mean2)**2.0_dp + (error1 / mean1)**2.0_dp &
- 2.0_dp * covariance / ((size(this)) * mean1 * mean2))
final_eie = final_error / sqrt(2.0_dp * (size(this) - 1))
write(iunit, *) size(that), mean_proje, final_error, final_eie
blocking_events = int(log(real(length, dp)) / log(2.0_dp) - 1.0_dp)
do i = 1, blocking_events
call reblock_data(this, blocklength)
call analyze_data(this, mean1, error1, eie1) !means & error in num
call reblock_data(that, blocklength)
call analyze_data(that, mean2, error2, eie2) !means & error in denom
! now have the correct means, but incorrect errors
covariance = calc_covariance(that, this)
mean_proje = mean1 / mean2
associate (sqrt_kernel => (error2 / mean2)**2.0_dp &
+ (error1 / mean1)**2.0_dp &
- 2.0_dp * covariance / ((size(this)) * mean1 * mean2) &
)
final_error = abs(mean_proje) * sqrt(merge(sqrt_kernel, 0._dp, sqrt_kernel > 0))
end associate
final_eie = final_error / sqrt(2.0_dp * (size(this) - 1))
write(iunit, *) size(that), mean_proje, final_error, final_eie
end do
close(iunit)
end subroutine print_proje_blocks