# Getting my method of bisection to output the correct value

• Fortran
I am trying to write a program that calculates the root of chi-square. I am not getting the correct answer and I honestly am at my wits end trying to figure it out. I know my simp_p() method is returning the correct value, but for some reason my root_chisq() method is not giving me the correct values. It's giving -A- value, but it isn't the right one. Can somebody help me out?

Fortran:
program main

!change V, B, and P in ROOT_CHISQ subroutine to see what answers

implicit none

real :: A, B, V, final
integer :: N

a=0
b=12.5
v=10
n=1000000

call simp_p(A,B,N,V,final)
print *, final

call root_chisq()

end program main

subroutine simp_p(A,B,N,V,final)

real :: A,B ! endpoints and number of divisions
real :: DX, integral, x, place
real :: V ! degrees freedom
real :: answer, final ! return from subroutine
integer :: N, i
integer(selected_int_kind(r=12)) :: i12
! calculate partition width

DX = (B-A)/FLOAT(N)
!print *, "DX = ", DX
!print *, "V = ", V

! initialize integral variable
integral = 0

integral = integral + 0
integral = integral + 0
else
end if

integral = integral + 0
integral = integral + 0
else
end if

!now compute values of remaining partitions

do i=1, N-1
!print *, "i equals = ", i
x = A+i*DX
!print *, "xval = ", x

cycle
cycle
end if

if (MOD(i,2).eq.0) THEN

else if (mod(i,2).ne.0) THEN
end if

!print *, integral
!print *,
!print *,

end do

!print *, 'final: ', integral
integral = integral * DX/3

! EVERYTHING BELOW HAS BEEN CHECKED, SOMETHING IS WRONG
! WITH THE INTEGRATION
place = place**(-1)
!print *, "place = ", place
final= place*integral

print *, 'final ANSWER = ', final

end subroutine simp_p

real, intent(in) :: v, x
real :: pt1, pt2

pt1 = x**((v/2)-1)
pt2 = exp(-x/2)
return
end subroutine equation

integer :: n, i
n = 1
difference = 1

do i=1, 56140
n=n+1
end do

end subroutine calc_gamma

subroutine root_chisq()

! assuming root lies somewhere in the interval 0,26
real :: p = 0.75
real :: v=10
real :: error = 0.01
real :: b=100.0 ! highest possible cutoff value when v=10
integer :: N=1000
logical :: check=.false.
!first check value of function at two endpoints

else
print *, "KEEP SEARCHING"
end if

check=.true.
else

print *, "KEEP SEARCHING"
end if

do

m=(a+b)/2

exit

exit

exit

else

b=m
print *, "replacing b"
print *, "m = ", (a+b)/2
else
a=m
print *, "replacing a"
print *, "m = ", (a+b)/2
end if
end if

if (abs(b-a).LT.e) THEN
print *, "final answerb = ", b
print *, "final answera = ", a

print *, (a+b)/2
print *, b
exit
end if

end do
return

end subroutine root_chisq

Mark44
Mentor
I am trying to write a program that calculates the root of chi-square. I am not getting the correct answer and I honestly am at my wits end trying to figure it out. I know my simp_p() method is returning the correct value, but for some reason my root_chisq() method is not giving me the correct values. It's giving -A- value, but it isn't the right one.
How far off is the value from root_chisq()? Is it just a little bit off, or is it way off? One thing I notice is that you're using real for the variable type, which is only four bytes for a floating point number. Possibly using real(8) or doubleprecision would give you better results.

You say you are having problems with your root_chisq() routine. As written, this routine has no parameters, and uses a bunch of "magic" numbers. In the first executable line, you have call simp_p(a,a,N,V, answer), in which you are apparently integrating something from a to a. That should produce a value of zero, and might be the source of your problem.

Much better style for root_chisq() would be to pass whatever parameters it needs, rather than rely on global variables. Your instructor will likely mark your program down for this.

Finally, whatever system you're using almost certainly has some powerful debugging capabilities. It's a real time-saver to use a debugger when a program is running but producing incorrect values. I would strongly encourage you to get up to speed with at least the rudiments of your debugger.

Minor thing: In simp_p, you have this: place = place**(-1). While not incorrect, it's simpler and most likely more efficient to write place = 1.0/place. Raising to a power is much less efficient than just doing division.
sarah_frtrn said:
Can somebody help me out?

Fortran:
program main

!change V, B, and P in ROOT_CHISQ subroutine to see what answers

implicit none

real :: A, B, V, final
integer :: N

a=0
b=12.5
v=10
n=1000000

call simp_p(A,B,N,V,final)
print *, final

call root_chisq()

end program main

subroutine simp_p(A,B,N,V,final)

real :: A,B ! endpoints and number of divisions
real :: DX, integral, x, place
real :: V ! degrees freedom
real :: answer, final ! return from subroutine
integer :: N, i
integer(selected_int_kind(r=12)) :: i12
! calculate partition width

DX = (B-A)/FLOAT(N)
!print *, "DX = ", DX
!print *, "V = ", V

! initialize integral variable
integral = 0

integral = integral + 0
integral = integral + 0
else
end if

integral = integral + 0
integral = integral + 0
else
end if

!now compute values of remaining partitions

do i=1, N-1
!print *, "i equals = ", i
x = A+i*DX
!print *, "xval = ", x

cycle
cycle
end if

if (MOD(i,2).eq.0) THEN

else if (mod(i,2).ne.0) THEN
end if

!print *, integral
!print *,
!print *,

end do

!print *, 'final: ', integral
integral = integral * DX/3

! EVERYTHING BELOW HAS BEEN CHECKED, SOMETHING IS WRONG
! WITH THE INTEGRATION
place = place**(-1)
!print *, "place = ", place
final= place*integral

print *, 'final ANSWER = ', final

end subroutine simp_p

real, intent(in) :: v, x
real :: pt1, pt2

pt1 = x**((v/2)-1)
pt2 = exp(-x/2)
return
end subroutine equation

integer :: n, i
n = 1
difference = 1

do i=1, 56140
n=n+1
end do

end subroutine calc_gamma

subroutine root_chisq()

! assuming root lies somewhere in the interval 0,26
real :: p = 0.75
real :: v=10
real :: error = 0.01
real :: b=100.0 ! highest possible cutoff value when v=10
integer :: N=1000
logical :: check=.false.
!first check value of function at two endpoints

else
print *, "KEEP SEARCHING"
end if

check=.true.
else

print *, "KEEP SEARCHING"
end if

do

m=(a+b)/2

exit

exit

exit

else

b=m
print *, "replacing b"
print *, "m = ", (a+b)/2
else
a=m
print *, "replacing a"
print *, "m = ", (a+b)/2
end if
end if

if (abs(b-a).LT.e) THEN
print *, "final answerb = ", b
print *, "final answera = ", a

print *, (a+b)/2
print *, b
exit
end if

end do
return

end subroutine root_chisq

Tom.G
I agree with @Mark44 that using any available debugger would be both informative and useful for your understanding.

For instance, the program has multiple instances of: if (answer.ne.answer) THEN, which seems rather pointless!

Mark44
Mentor
For instance, the program has multiple instances of: if (answer.ne.answer) THEN, which seems rather pointless!
Good catch -- I didn't notice this, but then again, I didn't pore through the program. Good example of "dead code"; code that will never execute, that should be removed.

And the rest of the code is
Fortran:
if (answer.ne.answer) THEN
integral = integral + 0[/quote]

I agree with @Mark44 that using any available debugger would be both informative and useful for your understanding.

For instance, the program has multiple instances of: if (answer.ne.answer) THEN, which seems rather pointless!
This is to check for a value of NaN. A value of NaN does not equal itself.

Thank you!

Tom.G
Good catch -- I didn't notice this, but then again, I didn't pore through the program. Good example of "dead code"; code that will never execute, that should be removed.

And the rest of the code is
Fortran:
if (answer.ne.answer) THEN
integral = integral + 0[/quote]

The only time it is executed if there is a value of NaN which was causing errors before.

Last edited by a moderator:
Tom.G
This is to check for a value of NaN. A value of NaN does not equal itself.
Ahh! Sneaky. I withdraw my comment... but will surely try to remember yours!

Cheers,
Tom

How far off is the value from root_chisq()? Is it just a little bit off, or is it way off? One thing I notice is that you're using real for the variable type, which is only four bytes for a floating point number. Possibly using real(8) or doubleprecision would give you better results.

You say you are having problems with your root_chisq() routine. As written, this routine has no parameters, and uses a bunch of "magic" numbers. In the first executable line, you have call simp_p(a,a,N,V, answer), in which you are apparently integrating something from a to a. That should produce a value of zero, and might be the source of your problem.

Much better style for root_chisq() would be to pass whatever parameters it needs, rather than rely on global variables. Your instructor will likely mark your program down for this.

Finally, whatever system you're using almost certainly has some powerful debugging capabilities. It's a real time-saver to use a debugger when a program is running but producing incorrect values. I would strongly encourage you to get up to speed with at least the rudiments of your debugger.

Minor thing: In simp_p, you have this: place = place**(-1). While not incorrect, it's simpler and most likely more efficient to write place = 1.0/place. Raising to a power is much less efficient than just doing division.

1. I will try this.
2. Yes it does return a value of 0, but then I subtract p from it which gives a negative p value. When doing method of bisection, you check the value of the function at each endpoint. My function subtracts p outside of the integral, so p doesn't get subtracted in the simp_p subroutine.
3. I am going to switch it so that I pass in parameters .. I was just using this for the time being until I have the code working correctly.
4. I've never used a debugger before. I will try to figure out how to run one and try it out. Thanks for all of your help!

How far off is the value from root_chisq()? Is it just a little bit off, or is it way off? One thing I notice is that you're using real for the variable type, which is only four bytes for a floating point number. Possibly using real(8) or doubleprecision would give you better results.

You say you are having problems with your root_chisq() routine. As written, this routine has no parameters, and uses a bunch of "magic" numbers. In the first executable line, you have call simp_p(a,a,N,V, answer), in which you are apparently integrating something from a to a. That should produce a value of zero, and might be the source of your problem.

Much better style for root_chisq() would be to pass whatever parameters it needs, rather than rely on global variables. Your instructor will likely mark your program down for this.

Finally, whatever system you're using almost certainly has some powerful debugging capabilities. It's a real time-saver to use a debugger when a program is running but producing incorrect values. I would strongly encourage you to get up to speed with at least the rudiments of your debugger.

Minor thing: In simp_p, you have this: place = place**(-1). While not incorrect, it's simpler and most likely more efficient to write place = 1.0/place. Raising to a power is much less efficient than just doing division.

Oh, and for the amount my answer is off...

SUPPOSED TO BE: 11.4, getting 11.8
supposed to be: 12.6, getting 12.5
supposed to be 10.2, getting 11.0

So I figured out what is happening. I have been giving it a cutoff point (the value I'm starting with for B....) and every time I'm getting roughly half of what B was at the beginning. I'm not too sure what is going wrong... but at least I figured out the problem.

Mark44
Mentor
2. Yes it does return a value of 0, but then I subtract p from it which gives a negative p value.
But what's the point of this? If you subtract p from 0, you get -p.
sarah_frtrn said:
When doing method of bisection, you check the value of the function at each endpoint.
This is what I was asking about, at the beginning of the executable statements in root_chisq():
Fortran:
call simp_p(a,a,N,V, answer)
answera = answer-p
In the call to simp_p(), it seems to me that answer gets set to 0.0, so you print 0.0, and then -p, and then set answera to -p,
I don't understand why you're doing this -- it would be simpler to just set answera to -p.

I have no idea why you're doing any of this, and there are no comments to help me understand it. At least based on the name, your simp_p() routine should calculate the integral of some function between x = a and x = b, and put the result in answer.
sarah_frtrn said:
My function subtracts p outside of the integral, so p doesn't get subtracted in the simp_p subroutine.
Again, I have no idea why you are doing this.

Okay I will explain. Yes, simp_p() routine calculates the integral of the equation found in the equation routine.

Basically, I first wrote the simp_p() routine because, given a cutoff value (obtained from a table on the percentile values for the chi-squared distribution), I am able to calculate the area under the curve ---- giving me the p-value.

This routine worked fine.

Then, I tried to reverse it, given a p-value along with the degrees of freedom, what is the cutoff value?
This is what the subroutine root_chisq() is trying to calculate. By setting the equation equal to 0 (((subtracting p))), and finding the roots (through the method of bisection), I will then have my cutoff value.

(So I can thereby run a hypothesis test on fortran knowing only my significance level and degrees of freedom.

Now, this would all be great if I could actually get the program to work.

Mark44
Mentor
Mark44 said:
And the rest of the code is
Fortran:
if (answer.ne.answer) THEN
integral = integral + 0
The only time it is executed if there is a value of NaN which was causing errors before.
Wouldn't it be better to find out why you're getting NAN in the first place, than to include this kludge?

IMO, you would be better off rewriting your root_chisq() routine, which I quote below. I have slightly revised your formatting, and have added line numbers to be able to identify places I'm commenting on.
Fortran:
subroutine root_chisq()

! assuming root lies somewhere in the interval 0,26
real :: p = 0.75
real :: v=10
real :: error = 0.01
real :: b=100.0 ! highest possible cutoff value when v=10
integer :: N=1000               !! Line 10
logical :: check=.false.
!first check value of function at two endpoints

print *, "ANSWER FOUND"               !! Line 20
else
print *, "KEEP SEARCHING"
end if

check=.true.                                      !! Line 30
else
print *, "KEEP SEARCHING"
end if

do
m=(a+b)/2

exit

exit

exit                                                     !! Line 60

else

b=m
print *, "replacing b"
print *, "m = ", (a+b)/2
else
a=m                                                    !! Line 70
print *, "replacing a"
print *, "m = ", (a+b)/2
end if
end if

if (abs(b-a).LT.e) THEN
print *, "final answerb = ", b
print *, "final answera = ", a
!! Line 80

print *, (a+b)/2
print *, b
exit
end if

end do
return                                                      !! Line 90

end subroutine root_chisq

1. What does v represent?
2. error set to .01 might be too large
3. Line 14 -- calling simp_p(a, a, ...) seems like a complete waste. It's like making the computer evaluate ##\int_a^a f(t) dt##, which is just zero (assuming f is an integrable function).
Fortran:
call simp_p(a,a,N,V, answer)
answera = answer-p
As mentioned before, the first line sets answer to 0.0, then prints it, then prints -.75 (p's hard-coded value), then sets answera to -.75.
This code could be (and should be) pruned way back.
4. Line 40 -- here you're calling simp_p(a, a, ...) again, which should again result in answer being set to 0.0. Why is it there?
5. Lines 28 and 50 and 54 and 58 -- In line 28, you're comparing |answerb| to .0001, but in lines 50 and 54, you're comparing |answerm| and |answerb| to .001. And similar in line 58. Be consistent. Pick one value and make it a parameter. Your variable error (set to .01) seems a likely candidate to be a parameter. You should not be using all of these "magic" numbers in your code.
6. Line 77 -- "if (abs(b-a).LT.e) THEN" -- what is e? This variable is not declared in root_chisq(), nor anywhere else that I'm able to find. It's almost certainly some garbage value.
7. You have lots of print statements, which is good, since you're not familiar with a debugger. Print statements are a primitive way of debugging, but definitely effective if that's all you have to work with.

Last edited:
jim mcnamara