1
$\begingroup$

I'm trying to implement a rescale/normalization as is described in here by mathematical expressions in page 10, and code in the top part of page 16. The notebook is also available here. The objective of the normalization is for the elements, created in a recursive way, of a vector to never go below machine-size, and to be able to still be able to use the normalized values to compute some variables of interest.

Here is my adapted (non)working example:

n = 999;
alpha = RandomReal[{0, 1}, {2, n}, WorkingPrecision -> 20];
    (*the totals below are just for another rescalement. 
They don't seem to be the problem, since even without them the program gets the same error.*)
    bmat = (#/Total[#] &) /@ 
       RandomReal[{0, 1}, {2, n}, WorkingPrecision -> 20];
    miota = ((#/Total[#] &) /@ 
     RandomReal[{0, 1}, {1, 2}, WorkingPrecision -> 20])[[1]];
mtransprob = (#/Total[#] &) /@ 
   RandomReal[{0, 1}, {2, 2}, WorkingPrecision -> 20];

listalpha = Range[n];
cf = Compile[{{y, _Real}},

   alpha[[1, 1]] = miota[[1]]*bmat[[1, 1]];
   alpha[[2, 1]] = miota[[2]]*bmat[[2, 1]];

   aux = Transpose[mtransprob];
   ct = 1;

   While[ct <= n - 1,

    alpha[[1 ;; 2, ct + 1]] = (aux.alpha[[1 ;; 2, ct]])*
      bmat[[1 ;; 2, ct + 1]];

    listalpha[[ct + 1]] = 1/(Total[alpha[[1 ;; 2, ct + 1]]]);
    (*The rescalement is done here. 
It's just an average, that should make the elements of the vector alpha with value between ]0,1[ 
and prevent the values of the alpha recursion from going below machine-size. However it doesn't seem to work*)
    alpha[[1 ;; 2, ct + 1]] = 
     alpha[[1 ;; 2, ct + 1]]*listalpha[[ct + 1]];

    ct = ct + 1;

    ], CompilationTarget -> "C"];

So, if I run cf[1], which will call the function cf and update the matrix alpha, I get the error:

Power::infy: "Infinite expression 1/0.*10^-3 encountered. "

I get the same error, even if I remove the compile command. I think it has to do with the way the initial variables are initiated. But I don't understand why, since I'm applying the same reasoning as this tutorial...

Edit1: It seems if at some steps of the while I put //N the error doesn't show up anymore...

$\endgroup$

Your Answer

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

Browse other questions tagged or ask your own question.