1
votes

I am using Mathematica 9 Student Edition and I am having problems with recursion limit errors. I have written a procedure that works and involves numerical integration. If I run this procedure by itself and plug in the values I want to test, then this procedure works just fine. The procedure is printed below. For anyone interested, it is supposed to transform two random variables, betaone and betazero into two other random variables gamma and rho and then numerically integrate a probability density function using these new random variables:

xtheta = .4; xx = .1;

c = 1/theta /. theta -> xtheta;

gamma = (Log[theta / (1 - theta)] - betazero - x*betaone) / 
   betaone;  
rho = Exp[ betazero + x*betaone] / (1 + 
    Exp[betazero + x*betaone]); 

JMatrix = {{D[gamma, betazero], 
   D[gamma, betaone]}, {D[rho, betazero], D[rho, betaone]}};

myJacobian = Det[JMatrix] ; gamma =.; rho =.; 

betazero = (1/
    gamma)*((x + gamma)*Log[rho/(1 - rho)] - x*Log[theta/(1 - theta)]);

betaone = (1/gamma)*(Log[theta/(1 - theta)] - Log[rho/(1 - rho)]);

finalPDF = 
  c * myJacobian  /. {betazero -> (1/
        gamma)*((x + gamma)*Log[rho/(1 - rho)] - 
        x*Log[theta/(1 - theta)]), 
    betaone -> (1/gamma)*(Log[theta/(1 - theta)] - 
        Log[rho/(1 - rho)])}; 

theta = xtheta; finalPDF2 = finalPDF /. {x -> xx};

n = NIntegrate[finalPDF2, {rho, 0, xtheta - .001},  {gamma, 0, 1}];

However, as soon as I add a Print[n] statement right after this procedure the code cannot run and I get the following error "$RecursionLimit::reclim: Recursion depth of 1024 exceeded". Can anyone explain why I get this error and how I can remedy it?

Additionally, I hope to use this procedure in a loop, so that I can run this procedure over and over again and change one of the variables slightly and output my results using the Print statement. I have tried using Do loops and For loops, but I get the same problem as before where I get a "$RecursionLimit::reclim: Recursion depth of 1024 exceeded" error message. Does anyone know if there is an advantage to using one loop over the other and if the problem with the looping is the same as my problem with the Print statement?

Thank you for your help!

1
youve assigned betazero/one then use them as patter substituions - agentp
Followup, the code runs fine on first pass, the error results when you re-run it, as betazero,betaone are now defined and you are taking derivatives with respct to those complicated expressions. Simply eliminating the assignments fixes it (or clear them at the top) - agentp
There is a dedicated StackExchange site for Mathematica now: mathematica.stackexchange.com -- please ask your future questions there. - Mr.Wizard

1 Answers

1
votes

For sake os providing an answer, here is a cleaned up version of your code. Note we avoid ever assigning to things we want to use as symbols for integration/derivatives.

f[xtheta_, xx_] :=
    Module[{theta, gamma, rho, betazero, betaone, x, sub, myJacobian, 
       JMatrix, finalPDF},
   sub = {
      gamma -> (Log[theta/(1 - theta)] - betazero - x*betaone)/betaone,
      rho -> Exp[betazero + x*betaone]/(1 + Exp[betazero + x*betaone])};
  JMatrix = {
    {D[gamma /. sub, betazero], D[gamma /. sub, betaone]},
    {D[rho /. sub, betazero],D[rho /. sub, betaone]}} ;
   myJacobian = Det[JMatrix];
   finalPDF = myJacobian/theta /.
     {theta -> xtheta,
     betazero -> (1/gamma)*((x + gamma)*Log[rho/(1 - rho)] - 
      x*Log[theta/(1 - theta)]), 
     betaone -> (1/gamma)*(Log[theta/(1 - theta)] - 
      Log[rho/(1 - rho)])};
   NIntegrate[finalPDF /. {x -> xx, theta -> xtheta},
      {rho, 0, xtheta - .001},
      {gamma, 0, 1}]];

n = f[.4, .1]

(* 0.247738 *)

The module construct makes this "safe" to loop, ie. guarantees no unintended side effects from one pass to the next:

Table[f[p, q], {p, {.2, .4, .6}}, {q, {.1, .2}}]

({{0.17888, 0.17888}, {0.247738, 0.247738}, {0.197697, 0.197697}})