Getting my method of bisection to output the correct value

[sarah_frtrn]

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?

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 mainsubroutine 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

call equation(A, V, answer)

!print *, "answer a = ", answer

if (answer.ne.answer) THEN
integral = integral + 0
else if (answer.ge.huge(i12)) THEN
integral = integral + 0
else
integral = integral + answer
end if

call equation(B, V, answer)
!print *, "answer b = ", answer

if (answer.ne.answer) THEN
integral = integral + 0
else if (answer.ge.huge(i12)) THEN
integral = integral + 0
else
integral = integral + answer
end if

!now compute values of remaining partitions

do i=1, N-1
!print *, "i equals = ", i
x = A+i*DX
!print *, "xval = ", x
call equation(x, v, answer )

if (answer.ne.answer) THEN
cycle
else if (answer.ge.huge(i12)) THEN
cycle
end if

if (MOD(i,2).eq.0) THEN
integral = integral + 2*answer

else if (mod(i,2).ne.0) THEN
INTEGRAL = INTEGRAL + 4*answer
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
call calc_gamma(v/2, answer)
place = (2**(v/2))*answer
place = place**(-1)
!print *, "place = ", place
final= place*integral

print *, 'final ANSWER = ', final
end subroutine simp_p
subroutine equation(x, v, answer)
real, intent(in) :: v, x
real, intent(out) :: answer
real :: pt1, pt2

pt1 = x**((v/2)-1)
pt2 = exp(-x/2)
answer = pt1*pt2
!print *, "answer from eq = ", answer
return
end subroutine equation
subroutine calc_gamma(x, answer)real :: x, answer, difference, previous_answer
integer :: n, i
answer = 1
n = 1
difference = 1
previous_answer = 1

do i=1, 56140
answer = answer * (((1/float(n)) + 1)**x)/((x/float(n))+ 1)
difference = abs(previous_answer - answer)
n=n+1
end doanswer = answer * (1/x)

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 :: a=0, answera, answerb, m, answerm
real :: b=100.0 ! highest possible cutoff value when v=10
real :: answer
integer :: N=1000
logical :: check=.false.
!first check value of function at two endpoints
call simp_p(a,a,N,V, answer)
print *, "answera = ", answer
print *, "answer - p = ", answer-p
answera = answer-p

if (abs(answera).LE.0.0001) THEN
  print *, "ANSWER FOUND"
else
  print *, "KEEP SEARCHING"
end if

call simp_p(a,b,n,v,answer)

answerb = answer-p
print *, "answer = ", answer
if (abs(answerb).LE.0.0001) THEN
  print *, "ANSWER FOUND"
  check=.true.
else

  print *, "KEEP SEARCHING"
end if

do

m=(a+b)/2
call simp_p(a,m,n,v,answer)
answerm = answer-p

call simp_p(a,a,n,v,answer)
answera = answer-p

call simp_p(a,b,n,v,answer)
answerb = answer-p

print *, "answer a = ", answera
print *, "answer b = ", answerb
print *, "answer m = ", answerm

if(abs(answerm).le.0.001) THEN
print *, "FINAL ANSWERM = ", answerm
exit

else if (abs(answerb).le.0.001) THEN
print *, "FINAL ANSWERB = ", answerb
exit

else if (abs(answera).le.0.001) THEN
print *, "FINAL ANSWERA = ", answera
exit

else

if ((answera.lt.0).AND.(answerm.GT.0)) THEN
  b=m
  answerb=answerm
  print *, "replacing b"
  print *, "m = ", (a+b)/2
else
a=m
answera=answerm
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 *, "answerb = ", answerb
print *, "answera = ", answera

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

end do
return

end subroutine root_chisq

 

[Mark44]

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.

Can somebody help me out?

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 mainsubroutine 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

call equation(A, V, answer)

!print *, "answer a = ", answer

if (answer.ne.answer) THEN
integral = integral + 0
else if (answer.ge.huge(i12)) THEN
integral = integral + 0
else
integral = integral + answer
end if

call equation(B, V, answer)
!print *, "answer b = ", answer

if (answer.ne.answer) THEN
integral = integral + 0
else if (answer.ge.huge(i12)) THEN
integral = integral + 0
else
integral = integral + answer
end if

!now compute values of remaining partitions

do i=1, N-1
!print *, "i equals = ", i
x = A+i*DX
!print *, "xval = ", x
call equation(x, v, answer )

if (answer.ne.answer) THEN
cycle
else if (answer.ge.huge(i12)) THEN
cycle
end if

if (MOD(i,2).eq.0) THEN
integral = integral + 2*answer

else if (mod(i,2).ne.0) THEN
INTEGRAL = INTEGRAL + 4*answer
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
call calc_gamma(v/2, answer)
place = (2**(v/2))*answer
place = place**(-1)
!print *, "place = ", place
final= place*integral

print *, 'final ANSWER = ', final
end subroutine simp_p
subroutine equation(x, v, answer)
real, intent(in) :: v, x
real, intent(out) :: answer
real :: pt1, pt2

pt1 = x**((v/2)-1)
pt2 = exp(-x/2)
answer = pt1*pt2
!print *, "answer from eq = ", answer
return
end subroutine equation
subroutine calc_gamma(x, answer)real :: x, answer, difference, previous_answer
integer :: n, i
answer = 1
n = 1
difference = 1
previous_answer = 1

do i=1, 56140
answer = answer * (((1/float(n)) + 1)**x)/((x/float(n))+ 1)
difference = abs(previous_answer - answer)
n=n+1
end doanswer = answer * (1/x)

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 :: a=0, answera, answerb, m, answerm
real :: b=100.0 ! highest possible cutoff value when v=10
real :: answer
integer :: N=1000
logical :: check=.false.
!first check value of function at two endpoints
call simp_p(a,a,N,V, answer)
print *, "answera = ", answer
print *, "answer - p = ", answer-p
answera = answer-p

if (abs(answera).LE.0.0001) THEN
  print *, "ANSWER FOUND"
else
  print *, "KEEP SEARCHING"
end if

call simp_p(a,b,n,v,answer)

answerb = answer-p
print *, "answer = ", answer
if (abs(answerb).LE.0.0001) THEN
  print *, "ANSWER FOUND"
  check=.true.
else

  print *, "KEEP SEARCHING"
end if

do

m=(a+b)/2
call simp_p(a,m,n,v,answer)
answerm = answer-p

call simp_p(a,a,n,v,answer)
answera = answer-p

call simp_p(a,b,n,v,answer)
answerb = answer-p

print *, "answer a = ", answera
print *, "answer b = ", answerb
print *, "answer m = ", answerm

if(abs(answerm).le.0.001) THEN
print *, "FINAL ANSWERM = ", answerm
exit

else if (abs(answerb).le.0.001) THEN
print *, "FINAL ANSWERB = ", answerb
exit

else if (abs(answera).le.0.001) THEN
print *, "FINAL ANSWERA = ", answera
exit

else

if ((answera.lt.0).AND.(answerm.GT.0)) THEN
  b=m
  answerb=answerm
  print *, "replacing b"
  print *, "m = ", (a+b)/2
else
a=m
answera=answerm
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 *, "answerb = ", answerb
print *, "answera = ", answera

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]

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

if (answer.ne.answer) THEN
integral = integral + 0[/quote]

 

[sarah_frtrn]

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!

 

[sarah_frtrn]

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

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.

 

[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

 

[sarah_frtrn]

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!

 

[sarah_frtrn]

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]

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.

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():

call simp_p(a,a,N,V, answer)
print *, "answera = ", answer
print *, "answer - p = ", answer-p
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.

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.

 

[sarah_frtrn]

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]

And the rest of the code is

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.

subroutine root_chisq()

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

call simp_p(a,a,N,V, answer)
print *, "answera = ", answer
print *, "answer - p = ", answer-p
answera = answer-p

if (abs(answera).LE.0.0001) THEN
   print *, "ANSWER FOUND"               ! Line 20
else
   print *, "KEEP SEARCHING"
end if

call simp_p(a,b,n,v,answer)
answerb = answer-p
print *, "answer = ", answer
if (abs(answerb).LE.0.0001) THEN
   print *, "ANSWER FOUND"
   check=.true.                                      ! Line 30
else
   print *, "KEEP SEARCHING"
end if

do
   m=(a+b)/2
   call simp_p(a,m,n,v,answer)
   answerm = answer-p

   call simp_p(a,a,n,v,answer)              ! Line 40
   answera = answer-p

   call simp_p(a,b,n,v,answer)
   answerb = answer-p

   print *, "answer a = ", answera
   print *, "answer b = ", answerb
   print *, "answer m = ", answerm

   if(abs(answerm).le.0.001) THEN         ! Line 50
       print *, "FINAL ANSWERM = ", answerm
       exit

  else if (abs(answerb).le.0.001) THEN
      print *, "FINAL ANSWERB = ", answerb
      exit

  else if (abs(answera).le.0.001) THEN
      print *, "FINAL ANSWERA = ", answera
      exit                                                     ! Line 60

  else

  if ((answera.lt.0).AND.(answerm.GT.0)) THEN
      b=m
      answerb=answerm
      print *, "replacing b"
      print *, "m = ", (a+b)/2
  else
     a=m                                                    ! Line 70
     answera=answerm
     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 *, "answerb = ", answerb
   print *, "answera = ", answera

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

end do
return                                                      ! Line 90

end subroutine root_chisq

Questions and comments
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).

call simp_p(a,a,N,V, answer)
print *, "answera = ", answer
print *, "answer - p = ", answer-p
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.