1
$\begingroup$

I am new to Mathematica and wanted to solve an equation that depends on parameters. The code is below

    theta = 0.3;
    v = 0.05;
    alpha = 1;

    q := InverseCDF[ GammaDistribution[ alpha, 1], (1 - theta)]
    rk[k_, r_] := r k
    rkplus1[k_, r_] := r k + 1

    Lfun[k_, r_] := 
      (1/rkplus1[k, r])^alpha Exp[ rk[k,r] q] CDF[ GammaDistribution[ alpha, 
                                                                     1/(rk[k, r] +1)],
                                                   q]

    Lplustheta[k_, r_] := theta + Lfun[k, r]
    logLplustheta[k_, r_] := Log[ Lplustheta[k, r]]
    LHS[k_, r_] := rk[k, r] q Lfun[k, r]/Lplustheta[k, r] - logLplustheta[k, r  ]
    RHS[r_] := r v
    r := 1;

   While[ r < 100,
         { While[ LHS[ k 0.00001, r 0.05] < RHS[r 0.05], a = k 0.00001; k++], 
           Print[ r 0.05, " ", v, " ", theta, " ",
                  Lfun[k 0.00001, r 0.05], " ", a, " ", 
                  a/Lplustheta[ k 0.00001, r 0.05], " ", 
                  a - (a/Lplustheta[k 0.00001, r 0.05])]; k := 1}; r++]

The code works but is extremely slow. Is there a faster way to accomplish what I trying to do?
Essentially solve RHS - LHS == 0 for k, for a range of r, and print?

$\endgroup$
  • $\begingroup$ To quick things. 1) Don't use SetDelayed ( := ) to compute q, use Set instead. q is a constant and doesn't need to evaluated over and over again as you iterate through your loop. Define it once and for all with q = InverseCDF[GammaDistribution[alpha, 1], (1 - theta)]. 2) Don't write slow procedural code to solve your equation. Use one of Mathematica's built-in solvers to do it for you. $\endgroup$ – m_goldberg Aug 14 '14 at 15:11
3
$\begingroup$

It's your inner while loop that is causing the trouble. I refactored your code a bit also. Basically your value of k decreases after the first run. So you only have to do this whole run (over k) once. This is much faster:

θ = 0.3;
v = 0.05;
α = 1;
r = 1;
k = 1;
q = InverseCDF[GammaDistribution[α, 1], 1 - θ];
Lfun[kr_] := 
 Exp[kr q]/(kr + 1)^α CDF[
   GammaDistribution[α, 1/(kr + 1)], q]
LHS[kr_] := 
 Block[{Z = Lfun@kr}, (kr*q*Z)/(θ + Z) - Log[θ + Z]]
While[r < 100, 
 If[1 == r, While[LHS[k*0.00001*r*0.05] < v*r*0.05, k++];, 
  While[LHS[k*0.00001*r*0.05] > v*r*0.05, --k];
  ++k;];
 a = (k - 1)*0.00001;
 σ = Lfun[k*0.00001*r*0.05];
 Print[r*0.05, " ", v, " ", θ, " ", σ, " ", a, " ", 
  a/(θ + σ), " ", a - (a/(θ + σ))];
 r++]
$\endgroup$
  • $\begingroup$ Thank you very much. This works beautifully. $\endgroup$ – Shu .Anu Aug 14 '14 at 16:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.