! ch13_bootstrap_fit.f95
! Companion code for Chapter 13: Optimization and Data Fitting.
!
! Pairs-bootstrap uncertainty estimates for the fit y = a0 + a1*x.
! Resample complete (x,y) pairs with replacement and refit each sample.
! Assumes independent sampling units. For fixed experimental x values with
! a specified noise model, consider a residual or parametric bootstrap.
! Correlated observations require a scheme preserving their dependence.
! Standard errors are not confidence intervals and exclude systematic errors
! and model inadequacy. Fits are unweighted. Console output only.
!
! Compile (Fortran 95):
! gfortran -std=f95 -Wall -Wextra -O2 ch13_bootstrap_fit.f95 -o ch13_bootstrap_fit
! Run: ./ch13_bootstrap_fit (Windows: ch13_bootstrap_fit.exe)

module bootstrap_methods
    implicit none
    private
    integer, parameter :: dp=selected_real_kind(15,307)
    public :: dp, bootstrap_line, initialize_random, normal_random
contains
    subroutine line_coefficients(x,y,coeff)
        real(dp), intent(in) :: x(:),y(:)
        real(dp), intent(out) :: coeff(2)
        real(dp) :: xbar,ybar,Sxx,dx(size(x))
        if (size(x)<3 .or. size(y)/=size(x)) stop 'Supply matching vectors with at least three points.'
        if (.not. all(abs(x)<=huge(1.0_dp)) .or. .not. all(abs(y)<=huge(1.0_dp))) stop 'Data must be finite.'
        xbar=sum(x)/real(size(x),dp)
        ybar=sum(y)/real(size(y),dp)
        dx=x-xbar
        Sxx=sum(dx**2)
        if (Sxx<=0.0_dp .or. .not. (Sxx<=huge(1.0_dp))) then
            write(*,'(a)') 'A data set has no usable x spread; its slope is undefined.'
            stop 'Pairs bootstrap may be unsuitable for small or repeated-x data sets.'
        end if
        coeff(2)=sum(dx*(y-ybar))/Sxx
        coeff(1)=ybar-coeff(2)*xbar
        if (.not. all(abs(coeff)<=huge(1.0_dp))) stop 'Fit overflow: rescale measurements.'
    end subroutine line_coefficients

    subroutine bootstrap_line(x,y,B,estimate,samples,standard_errors,bias,covariance)
        real(dp), intent(in) :: x(:),y(:)
        integer, intent(in) :: B
        real(dp), intent(out) :: estimate(2),samples(:,:),standard_errors(2),bias(2),covariance(2,2)
        real(dp) :: xb(size(x)),yb(size(x)),u,means(2),delta(2)
        integer :: n,rep,i,j,k
        n=size(x)
        if (B<2 .or. size(samples,1)/=2 .or. size(samples,2)/=B) stop 'Invalid bootstrap sample array.'
        call line_coefficients(x,y,estimate)
        do rep=1,B
            do i=1,n
                call random_number(u)
                j=1+int(n*u)
                ! One index selects both coordinates of the observed pair.
                xb(i)=x(j)
                yb(i)=y(j)
            end do
            call line_coefficients(xb,yb,samples(:,rep))
        end do
        means=sum(samples,dim=2)/real(B,dp)
        bias=means-estimate
        covariance=0.0_dp
        do rep=1,B
            delta=samples(:,rep)-means
            do i=1,2
                do k=1,2
                    covariance(i,k)=covariance(i,k)+delta(i)*delta(k)
                end do
            end do
        end do
        covariance=covariance/real(B-1,dp)
        standard_errors(1)=sqrt(covariance(1,1))
        standard_errors(2)=sqrt(covariance(2,2))
    end subroutine bootstrap_line
    subroutine initialize_random(seed)
        integer, intent(in) :: seed
        integer, allocatable :: seeds(:)
        integer :: n,i
        call random_seed(size=n)
        allocate(seeds(n))
        do i=1,n
            seeds(i)=1+modulo(modulo(seed,1000000)+37*modulo(i,1000000),1000000)
        end do
        call random_seed(put=seeds)
        deallocate(seeds)
        ! Repeatable on the same compiler/runtime; sequences depend on compiler.
    end subroutine initialize_random

    function normal_random() result(z)
        real(dp) :: z,u1,u2
        do
            call random_number(u1)
            if (u1>0.0_dp) exit
        end do
        call random_number(u2)
        z=sqrt(-2.0_dp*log(u1))*cos(2.0_dp*acos(-1.0_dp)*u2)
    end function normal_random

end module bootstrap_methods

program ch13_bootstrap_fit
    use bootstrap_methods, only: dp, bootstrap_line, initialize_random, normal_random
    implicit none
    integer, parameter :: n=25,B=2000,seed=1
    real(dp) :: x(n),y(n),estimate(2),samples(2,B),standard_errors(2),bias(2),covariance(2,2)
    integer :: i
    ! integer :: clock_seed  ! Uncomment for clock-based seeding.

    call initialize_random(seed)
    ! Alternative: comment out the line above, uncomment the clock_seed
    ! declaration above, and uncomment both lines below.
    ! call system_clock(count=clock_seed)
    ! call initialize_random(clock_seed)

    ! Example data: intercept 2, slope 0.5, noise standard deviation 0.15.
    do i=1,n
        x(i)=real(i-1,dp)/real(n-1,dp)
        y(i)=2.0_dp+0.5_dp*x(i)+0.15_dp*normal_random()
    end do
    ! Replace the loop above with measured x,y data. For a plain two-column
    ! file containing exactly n rows and no header, use:
    ! open(unit=10,file='measurements.txt',status='old',action='read')
    ! do i=1,n
    !     read(10,*) x(i),y(i)
    ! end do
    ! close(10)

    call bootstrap_line(x,y,B,estimate,samples,standard_errors,bias,covariance)
    write(*,'(/,a)') 'Pairs bootstrap for y = a0 + a1*x'
    write(*,'(a,i0)') 'Measurements:      ',n
    write(*,'(a,i0)') 'Bootstrap samples: ',B
    write(*,'(a)') 'Parameter order: intercept, slope'
    write(*,'(a,2es20.12)') 'Original estimates: ',estimate
    write(*,'(a,2es20.12)') 'Bootstrap SE:       ',standard_errors
    write(*,'(a,2es20.12)') 'Estimated bias:     ',bias
    write(*,'(a)') 'Bootstrap covariance matrix:'
    do i=1,2
        write(*,'(2es20.12)') covariance(i,:)
    end do
    write(*,'(a)') 'SE is the sample standard deviation of bootstrap estimates.'
    write(*,'(a)') 'Bias is the bootstrap mean minus the original estimate.'
    write(*,'(a)') 'Increase B to check stability of the uncertainty estimates.'
end program ch13_bootstrap_fit