+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: HP Prime (/forum-5.html)
+--- Thread: longfloat library HP PRIME (/thread-17377.html)
longfloat library HP PRIME - robmio - 08-21-2021 03:43 PM
Hello everyone, after reading the following article: “Distribution of the largest root of a matrix for Roy’s test in multivariate analysis of variance”, written by Marco Chiani in the “Journal of Multivariate Analysis 143 (2016) 467–471”, I translated the algorithm on page 469 for HP PRIME, to obtain the “p” level of significance used in the multivariate analysis of variance according to Roy.
The algorithm works well on the HP PRIME, but, if I use high input values, the calculation becomes cumbersome, and the results may be wrong or undefined. The algorithm uses the incomplete Beta function and the Gamma function.
The results could be correct even with high input values if the HP PRIME possessed the “long-float library” of Giac-xCas.
I have read some “posts” on this topic, but I ask you this question again: will the HP PRIME be provided with the Giac-XCas “log-float library” in the future? Is it legitimate to hope?
Sincerely, Roberto Mioni
RE: longfloat library HP PRIME - robmio - 08-21-2021 05:48 PM
This is the algorithm:
Code:
#cas
CDF_Roy(s,m,n,theta):=
BEGIN
LOCAL A, ii, j, b, a, adzero, aa, cc;
A:=MAKEMAT(0,s,s);
b:=MAKELIST(0,x,1,s);
cc(x):=sqrt(π)*Gamma((2*m+2*n+s+x+2)/
2)/(Gamma((2*m+x+1)/2)*Gamma((2*n+x+
1)/2)*Gamma(x/2));
FOR ii FROM 1 TO s DO
b:=REPLACE(b,ii,{(Beta(m+ii,n+1,theta)
^2)/2});
FOR j FROM ii TO s-1 DO
b:=REPLACE(b,j+1,{(m+j)/(m+j+n+1)*
b(j)-Beta(2*m+ii+j,2*n+2,theta)/
(m+j+n+1)});
a:=(Beta(m+ii,n+1,theta)*Beta(m+j+1,
n+1,theta)-2*b(j+1))*cc(ii)*
cc(j+1);
A:=REPLACE(A,{ii,j+1},[[a]]);
END;
END;
aa:={};
FOR ii FROM 1 TO s DO
aa:=CONCAT(aa,{Beta(m+ii,n+1,theta)*
cc(ii)});
END;
aa:=ListToMat(aa);
adzero:=MAKELIST(0,x,1,s+1,1);
adzero:=ListToMat(adzero);
IF odd(s)==1 THEN
A:=ADDCOL(A,aa,s+1);
A:=ADDROW(A,adzero,s+1);
END;
A:=A-TRN(A);
RETURN √(DET(A));
END;
#endRE: longfloat library HP PRIME - robmio - 08-23-2021 04:17 AM
Thanks for the reply. Part of the problem is related to the high values to use with the Gamma function. One method is to use the LogGamma function (see the code below) and translate the result with Exp (LogGamm). It works. There remains the problem of the incomplete Beta function, which cannot be "transformed" into other functions (at least I think).
Best regards, Roberto.
Code:
#cas
LogGamma(ασ):=
BEGIN
RETURN int(Psi(x),x,1,ασ);
END;
#end