For a statistics application, I need to evaluate tiny chi-square probabilities. But when this becomes smaller than matlab's double-precision limit realmin (around 1e-308), it returns 0.
x=50;
p=chi2cdf(x^2,4,'upper')
>> 0
All I need is to store the order of magnitude of the number and a few leading digits. I know that using vpa you can store tiny numbers smaller than realmin to arbitrary precision:
y = vpa('1e-5000')
>> 1.0e-5000
log10(y)
>> -5000.0
But this doesn't seem to work with chi2cdf:
x = vpa('50');
chi2cdf(x^2,4,'upper')
throws an error.
Does this mean that some matlab functions, such as chi2cdf, were simply not designed to work with variable-precision inputs? Is there any other way around?
It looks like chi2cdf does not support symbolic inputs. The chi-squared distribution function can be expressed in terms of the regularized incomplete gamma function, which is gammainc in Matlab, but that one does not support symbolic inputs either.
So you may have to define the complementary distribution function manually as a symbolic integral by means of int, using only functions that accept symbolic inputs:
x = 50; % input value
k = 4; % input value
d = 30; % number of digits
syms ks ts
result = int(ts^(ks/2-1)*exp(-ts/2)/2^(ks/2)/gamma(ks/2), ts, inf);
result = vpa(subs(subs(result, ts, x^2), ks, k), d)
gives
result =
1.69494234796096964467999437076e-540