Am 29.11.2015 um 23:25 schrieb Allin Cottrell
<cottrell(a)wfu.edu>:
> On Sun, 29 Nov 2015, Pindar Os wrote:
>
> Hi there, I’m still working on the restricted OLS regressions and in the meantime
even uploaded an early version of a function package that uses a SQP of GNU R and a
bootstrap routine to get coeffs and standard errors. In order to produce a nice show
case of the optimization problem I generate all combination of coefficients with a step
size of n. If there are 3 coeffs that are >=0 and should sum up to 1 the matrix can be
easily generated via loop. Since gretl has a special "for loop" one could use
this loop version in order to fill the percentages in the matrix in one step. However,
there is a bug. Have a look. I’m using the current x64 snapshot for windows.
> Cheers and a nice first Sunday in Advent Leon
>
> <hansl>
> scalar intN_1_a = 50 # step size of 2 percent
> scalar intStep = 1/intN_1_a
> eval intStep
> matrix mCombis = NA
> matrix mCombisALL = NA
>
> set stopwatch
> loop for i=0..intN_1_a --quiet
> loop for j=0..intN_1_a --quiet
> if i+j<=intN_1_a
> if i+j == 0
> mCombisALL = i ~ j ~ (intN_1_a-i-j)
> else
> mCombisALL |= i ~ j ~ (intN_1_a-i-j)
> endif
> endif
> endloop
> endloop
> eval $stopwatch/60
> eval rows(mCombisALL)
> # correct number of rows since
> scalar intSum3 = (intN_1_a)*3 + (1+intN_1_a-3+1)*(1+intN_1_a-3)/2
> # for three coeefs I found this closed formula, with more it’s another…
> eval intSum3
>
>
> set stopwatch
> loop for (i=0; i<=1; i+=intStep) --quiet
> loop for (j=0; j<=1; j+=intStep) --quiet
> if i+j<=1
> if i+j == 0
> mCombis = i ~ j ~ (1-i-j)
> else
> mCombis |= i ~ j ~ (1-i-j)
> endif
> endif
> endloop
> endloop
> eval rows(mCombis)
> eval $stopwatch
> # correct number
> scalar intSum3 = (intN_1_a)*3 + (1+intN_1_a-3+1)*(1+intN_1_a-3)/2
> eval intSum3
> # wrong, WHY? The reason has to do with the.xx999 values and the missing 1s.
>
> <hansl>
It's discussed in the chapter on loops in the User's Guide. You need to allow for
the fact that when your i and j (in the second example) "ought to" equal 1 if
calculation were done with infinite precision, they're likely to be slightly off.
Hence you want something like
loop for (i=0; i<1.001; i+=intStep)
loop for (j=0; j<1.001; j+=intStep)
if i+j < 1.001
...
Ok, thanks for the hint!
Best
Leon
For reference, at the gretl console:
? printf "%.18g\n", 50*(1/50)
1.00000000000000022
Allin Cottrell
_______________________________________________
Gretl-users mailing list
Gretl-users(a)lists.wfu.edu
http://lists.wfu.edu/mailman/listinfo/gretl-users