diff --git a/Generalized chi-square distribution.mltbx b/Generalized chi-square distribution.mltbx new file mode 100644 index 0000000..ba427a9 Binary files /dev/null and b/Generalized chi-square distribution.mltbx differ diff --git a/demos.xml b/demos.xml new file mode 100644 index 0000000..aa16346 --- /dev/null +++ b/demos.xml @@ -0,0 +1,17 @@ + + + Generalized chi-square distribution + toolbox + HelpIcon.DEMOS + + The generalized chi-square variable is a quadratic form of a normal variable, or equivalently, a linear sum of independent non-central chi-square variables and a normal variable. Try the Getting Started guide for a quick demo of all the functions. + + + + + other + GettingStarted + doc/html/GettingStarted.html + + + \ No newline at end of file diff --git a/doc/GettingStarted.html b/doc/GettingStarted.html new file mode 100644 index 0000000..1d65d31 --- /dev/null +++ b/doc/GettingStarted.html @@ -0,0 +1,109 @@ + +Generalized chi-square distribution · Getting started

Generalized chi-square distribution · Getting started

The generalized chi-square variable is a quadratic form of a normal variable, or equivalently, a linear sum of independent non-central chi-square variables and a normal variable.
Look into each function code or type help functionname for more features and documentation.

Calculate mean and variance

% gx2 parameters
lambda=[1 -10 2];
m=[1 2 3];
delta=[2 3 7];
sigma=5;
c=10;
[mu,v]=gx2stat(lambda,m,delta,sigma,c)
mu = -17
v = 1771

Generate random samples

r=gx2rnd(lambda,m,delta,sigma,c,[1 1e4]);

Calculate pdf and cdf

x=[10 25];
f=gx2pdf(x,lambda,m,delta,sigma,c)
f = 1×2
0.0121 0.0088 +
p=gx2cdf(x,lambda,m,delta,sigma,c,'AbsTol',0,'RelTol',1e-4)
p = 1×2
0.7150 0.8790 +

Plot pdf and sample histogram

fplot(@(x) gx2pdf(x,lambda,m,delta,sigma,c))
hold on; histogram(r,'normalization','pdf','displaystyle','stairs')
xline(mu,'-',{'\mu \pm \sigma'},'labelorientation','horizontal');
xline(mu-sqrt(v),'-'); xline(mu+sqrt(v),'-');
xlim([mu-3*sqrt(v),mu+3*sqrt(v)]); ylim([0 .015]); ylabel 'pdf'

Plot cdf

figure; fplot(@(x) gx2cdf(x,lambda,m,delta,sigma,c));
xline(mu,'-',{'\mu \pm \sigma'},'labelorientation','horizontal');
xline(mu-sqrt(v),'-'); xline(mu+sqrt(v),'-');
xlim([mu-3*sqrt(v),mu+3*sqrt(v)]); ylim([0 1]); xlabel x; ylabel 'cdf'

Calculate inverse cdf

x=gx2inv([0.5 0.9],lambda,m,delta,sigma,c)
x = 1×2
-8.7657 27.5320 +

Distribution of quadratic form of a normal variable

% normal parameters
mu=[1;2]; % mean
v=[2 1; 1 3]; % covariance matrix
q(\mathbf{x})=(x_1+x_2)^2-x_1-1 = [x1;x2]'*[1 1; 1 1]*[x1;x2] + [-1;0]'*[x1;x2] -1
quad.q2=[1 1; 1 1];
quad.q1=[-1;0];
quad.q0=-1;
% get gx2 parameters corr. to this quadratic form
[lambda,m,delta,sigma,c]=gx2_params_norm_quad(mu,v,quad)
lambda = 7.0000
m = 1
delta = 1.1086
sigma = 0.8452
c = -0.7602
p(q(x)<3)
p=gx2cdf(3,lambda,m,delta,sigma,c)
p = 0.3351
+
+ \ No newline at end of file diff --git a/doc/GettingStarted.mlx b/doc/GettingStarted.mlx new file mode 100644 index 0000000..d0496df Binary files /dev/null and b/doc/GettingStarted.mlx differ diff --git a/doc/html/GettingStarted.html b/doc/html/GettingStarted.html new file mode 100644 index 0000000..06a576a --- /dev/null +++ b/doc/html/GettingStarted.html @@ -0,0 +1,124 @@ + +Generalized chi-square distribution · Getting started

Generalized chi-square distribution · Getting started

The generalized chi-square variable is a quadratic form of a normal variable, or equivalently, a linear sum of independent non-central chi-square variables and a normal variable.
For more features and documentation, type:
doc gx2stat
doc gx2rnd
doc gx2cdf
doc gx2cdf_davies
doc gx2cdf_imhof
doc gx2cdf_ruben
doc gx2pdf
doc gx2inv
doc gx2_params_norm_quad
If you use this toolbox, please cite: A method to integrate and classify normal distributions.
Bugs/questions/comments to abhranil.das@utexas.edu.
Abhranil Das, Center for Perceptual Systems, University of Texas at Austin

Calculate mean and variance

% gx2 parameters
lambda=[1 -10 2];
m=[1 2 3];
delta=[2 3 7];
sigma=5;
c=10;
[mu,v]=gx2stat(lambda,m,delta,sigma,c)
mu = -17
v = 1771

Generate random samples

r=gx2rnd(lambda,m,delta,sigma,c,[1 1e4]);

Calculate pdf and cdf

x=[10 25];
f=gx2pdf(x,lambda,m,delta,sigma,c)
f = 1×2
0.0121 0.0088 +
p=gx2cdf(x,lambda,m,delta,sigma,c,'AbsTol',0,'RelTol',1e-4)
p = 1×2
0.7150 0.8790 +

Plot pdf and sample histogram

fplot(@(x) gx2pdf(x,lambda,m,delta,sigma,c))
hold on; histogram(r,'normalization','pdf','displaystyle','stairs')
xline(mu,'-',{'\mu \pm \sigma'},'labelorientation','horizontal');
xline(mu-sqrt(v),'-'); xline(mu+sqrt(v),'-');
xlim([mu-3*sqrt(v),mu+3*sqrt(v)]); ylim([0 .015]); ylabel 'pdf'

Plot cdf

figure; fplot(@(x) gx2cdf(x,lambda,m,delta,sigma,c));
xline(mu,'-',{'\mu \pm \sigma'},'labelorientation','horizontal');
xline(mu-sqrt(v),'-'); xline(mu+sqrt(v),'-');
xlim([mu-3*sqrt(v),mu+3*sqrt(v)]); ylim([0 1]); xlabel x; ylabel 'cdf'

Calculate inverse cdf

x=gx2inv([0.5 0.9],lambda,m,delta,sigma,c)
x = 1×2
-8.7657 27.5320 +

Distribution of quadratic form of a normal variable

% normal parameters
mu=[1;2]; % mean
v=[2 1; 1 3]; % covariance matrix
q(\mathbf{x})=(x_1+x_2)^2-x_1-1 = [x1;x2]'*[1 1; 1 1]*[x1;x2] + [-1;0]'*[x1;x2] -1
quad.q2=[1 1; 1 1];
quad.q1=[-1;0];
quad.q0=-1;
% get gx2 parameters corr. to this quadratic form
[lambda,m,delta,sigma,c]=gx2_params_norm_quad(mu,v,quad)
lambda = 7.0000
m = 1
delta = 1.1086
sigma = 0.8452
c = -0.7602
p(q(x)<3)
p=gx2cdf(3,lambda,m,delta,sigma,c)
p = 0.3351
+
+ \ No newline at end of file diff --git a/gx2.prj b/gx2.prj new file mode 100644 index 0000000..dccd162 --- /dev/null +++ b/gx2.prj @@ -0,0 +1,152 @@ + + + Generalized chi-square distribution + Abhranil Das + abhranil.das@utexas.edu + The University of Texas at Austin + Compute the statistics, pdf, cdf, inverse cdf and random numbers of the generalized chi-square distribution. + The generalized chi-square variable is a quadratic form of a normal variable, or equivalently, a linear sum of independent non-central chi-square variables and a normal variable. Try the Getting Started guide for a quick demo of all the functions. + ${PROJECT_ROOT}\gx2_icon.png + 1.6.1 + ${PROJECT_ROOT}\Generalized chi-square distribution.mltbx + + Statistics and Machine Learning Toolbox + Symbolic Math Toolbox + + + 19 + 15 + + + 12.0 + 8.6 + + + 6f9f2474-f422-495e-a245-dd928544dd14 + % List files contained in your toolbox folder that you would like to exclude +% from packaging. Excludes should be listed relative to the toolbox folder. +% Some examples of how to specify excludes are provided below: +% +% A single file in the toolbox folder: +% .svn +% +% A single file in a subfolder of the toolbox folder: +% example/.svn +% +% All files in a subfolder of the toolbox folder: +% example/* +% +% All files of a certain name in all subfolders of the toolbox folder: +% **/.svn +% +% All files matching a pattern in all subfolders of the toolbox folder: +% **/*.bak +% +gx2_icon.png + true + <?xml version="1.0" encoding="utf-8"?> +<examples> + <exampleCategory name="doc"> + <example name="GettingStarted" type="html"> + <file type="source">/doc/html/GettingStarted.html</file> + <file type="main">/doc/GettingStarted.mlx</file> + <file type="thumbnail"/> + </example> + </exampleCategory> +</examples> + + + + + ${PROJECT_ROOT}\doc\GettingStarted.mlx + + + true + + + + + + false + true + true + true + true + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ${PROJECT_ROOT} + + + ${PROJECT_ROOT}\demos.xml + ${PROJECT_ROOT}\doc + ${PROJECT_ROOT}\gx2_params_norm_quad.m + ${PROJECT_ROOT}\gx2cdf.m + ${PROJECT_ROOT}\gx2cdf_davies.m + ${PROJECT_ROOT}\gx2cdf_imhof.m + ${PROJECT_ROOT}\gx2cdf_ruben.m + ${PROJECT_ROOT}\gx2inv.m + ${PROJECT_ROOT}\gx2pdf.m + ${PROJECT_ROOT}\gx2rnd.m + ${PROJECT_ROOT}\gx2stat.m + + + + + + C:\Users\abhra\Google Drive\Geisler Lab\IntClassNorm\gx2\Generalized chi-square distribution.mltbx + + + + C:\Program Files\MATLAB\R2020b + + + + false + false + true + false + false + false + false + false + 10.0 + false + true + win64 + true + + + \ No newline at end of file diff --git a/gx2_icon.png b/gx2_icon.png new file mode 100644 index 0000000..63dd726 Binary files /dev/null and b/gx2_icon.png differ diff --git a/gx2_params_norm_quad.m b/gx2_params_norm_quad.m new file mode 100644 index 0000000..ea7710e --- /dev/null +++ b/gx2_params_norm_quad.m @@ -0,0 +1,61 @@ +function [lambda,m,delta,sigma,c]=gx2_params_norm_quad(mu,v,quad) + + % GX2_PARAMS_NORM_QUAD A quadratic form of a normal variable is distributed + % as a generalized chi-squared. This function takes the normal parameters + % and the quadratic coeffs and returns the parameters of the generalized + % chi-squared. + % + % Abhranil Das + % Center for Perceptual Systems, University of Texas at Austin + % If you use this code, please cite: + % A method to integrate and classify normal distributions. + % + % Example: + % mu=[1;2]; % mean + % v=[2 1; 1 3]; % covariance matrix + % % Say q(x)=(x1+x2)^2-x1-1 = [x1;x2]'*[1 1; 1 1]*[x1;x2] + [-1;0]'*[x1;x2] - 1: + % quad.q2=[1 1; 1 1]; + % quad.q1=[-1;0]; + % quad.q0=-1; + % + % [lambda,m,delta,sigma,c]=gx2_params_norm_quad(mu,v,quad) + % + % Required inputs: + % mu column vector of normal mean + % v normal covariance matrix + % quad struct with quadratic form coefficients: + % q2 matrix of quadratic coefficients + % q1 column vector of linear coefficients + % q0 scalar constant + % + % Outputs: + % lambda row vector of coefficients of the non-central chi-squares + % m row vector of degrees of freedom of the non-central chi-squares + % delta row vector of non-centrality paramaters (sum of squares of + % means) of the non-central chi-squares + % sigma sd of normal term + % c constant term + % + % See also: + % Interactive demos + % gx2rnd + % gx2stat + % gx2cdf + % gx2pdf + + % standardize the space + q2s=0.5*(quad.q2+quad.q2'); % symmetrize q2 + q2=sqrtm(v)*q2s*sqrtm(v); + q1=sqrtm(v)*(2*q2s*mu+quad.q1); + q0=mu'*q2s*mu+quad.q1'*mu+quad.q0; + + [R,D]=eig(q2); + d=diag(D)'; + b=(R'*q1)'; + + [lambda,~,ic]=unique(nonzeros(d)'); % unique non-zero eigenvalues + m=accumarray(ic,1)'; % total dof of each eigenvalue + delta=arrayfun(@(x) sum((b(d==x)).^2),lambda)./(4*lambda.^2); % total non-centrality for each eigenvalue + sigma=norm(b(~d)); + c=q0-sum(lambda.*delta); diff --git a/gx2cdf.m b/gx2cdf.m new file mode 100644 index 0000000..0e34cfe --- /dev/null +++ b/gx2cdf.m @@ -0,0 +1,78 @@ +function p=gx2cdf(x,lambda,m,delta,sigma,c,varargin) + + % GX2CDF Returns the cdf of a generalized chi-squared (a weighted sum of + % non-central chi-squares and a normal), using Ruben's [1962] method, + % Davies' [1973] method, or the native ncx2cdf, depending on the input. + % + % Abhranil Das + % Center for Perceptual Systems, University of Texas at Austin + % If you use this code, please cite: + % A method to integrate and classify normal distributions. + % + % Usage: + % p=gx2cdf(x,lambda,m,delta,sigma,c) + % p=gx2cdf(x,lambda,m,delta,sigma,c,'upper') + % p=gx2cdf(x,lambda,m,delta,sigma,c,'AbsTol',0,'RelTol',1e-7) + % + % Example: + % f=gx2pdf(25,[1 -5 2],[1 2 3],[2 3 7],5,0) + % + % Required inputs: + % x points at which to evaluate the CDF + % lambda row vector of coefficients of the non-central chi-squares + % m row vector of degrees of freedom of the non-central chi-squares + % delta row vector of non-centrality paramaters (sum of squares of + % means) of the non-central chi-squares + % sigma sd of normal term + % c constant term + % + % Optional positional input: + % 'upper' more accurate estimate of the complementary CDF when it's small + % + % Optional name-value inputs: + % 'AbsTol' absolute error tolerance for the output + % 'RelTol' relative error tolerance for the output + % The absolute OR the relative tolerance is satisfied. + % + % Output: + % p computed cdf + % + % See also: + % Interactive demos + % gx2cdf_davies + % gx2cdf_imhof + % gx2cdf_ruben + % gx2pdf + + parser = inputParser; + parser.KeepUnmatched = true; + addRequired(parser,'x',@(x) isreal(x)); + addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); + addRequired(parser,'m',@(x) isreal(x) && isrow(x)); + addRequired(parser,'delta',@(x) isreal(x) && isrow(x)); + addRequired(parser,'sigma',@(x) isreal(x) && isscalar(x)); + addRequired(parser,'c',@(x) isreal(x) && isscalar(x)); + addOptional(parser,'side','lower',@(x) strcmpi(x,'lower') || strcmpi(x,'upper') ); + addParameter(parser,'AbsTol',1e-10,@(x) isreal(x) && isscalar(x) && (x>=0)); + addParameter(parser,'RelTol',1e-6,@(x) isreal(x) && isscalar(x) && (x>=0)); + + parse(parser,x,lambda,m,delta,sigma,c,varargin{:}); + side=parser.Results.side; + + if ~sigma && length(unique(lambda))==1 + % native ncx2 fallback + if (sign(unique(lambda))==1 && strcmpi(side,'lower')) || (sign(unique(lambda))==-1 && strcmpi(side,'upper')) + p=ncx2cdf((x-c)/unique(lambda),sum(m),sum(delta)); + else + p=ncx2cdf((x-c)/unique(lambda),sum(m),sum(delta),'upper'); + end + elseif ~sigma && (all(lambda>0)||all(lambda<0)) + try + p=gx2cdf_ruben(x,lambda,m,delta,c,varargin{:}); + catch + p=gx2cdf_davies(x,lambda,m,delta,sigma,c,varargin{:}); + end + else + p=gx2cdf_davies(x,lambda,m,delta,sigma,c,varargin{:}); + end \ No newline at end of file diff --git a/gx2cdf_davies.m b/gx2cdf_davies.m new file mode 100644 index 0000000..c5bedfa --- /dev/null +++ b/gx2cdf_davies.m @@ -0,0 +1,98 @@ +function [p,flag]=gx2cdf_davies(x,lambda,m,delta,sigma,c,varargin) + + % GX2CDF_DAVIES Returns the cdf of a generalized chi-squared (a weighted + % sum of non-central chi-squares and a normal), using Davies' [1973] + % method. + % + % Abhranil Das + % Center for Perceptual Systems, University of Texas at Austin + % If you use this code, please cite: + % A method to integrate and classify normal distributions. + % + % Usage: + % p=gx2cdf_davies(x,lambda,m,delta,sigma,c) + % p=gx2cdf_davies(x,lambda,m,delta,sigma,c,'upper') + % p=gx2cdf_davies(x,lambda,m,delta,sigma,c,'AbsTol',0,'RelTol',1e-7) + % + % Example: + % p=gx2cdf_davies(25,[1 -5 2],[1 2 3],[2 3 7],5,0) + % + % Required inputs: + % x points at which to evaluate the cdf + % lambda row vector of coefficients of the non-central chi-squares + % m row vector of degrees of freedom of the non-central chi-squares + % delta row vector of non-centrality paramaters (sum of squares of + % means) of the non-central chi-squares + % sigma sd of normal term + % c constant term + % + % Optional positional input: + % 'upper' more accurate estimate of the complementary CDF when it's small + % + % Optional name-value inputs: + % 'AbsTol' absolute error tolerance for the output + % 'RelTol' relative error tolerance for the output + % The absolute OR the relative tolerance is satisfied. + % + % Outputs: + % p computed cdf + % flag =true if output was too close to 0 or 1 to compute exactly with + % default settings. Try stricter tolerances. + % + % See also: + % Interactive demos + % gx2cdf_imhof + % gx2cdf_ruben + % gx2cdf + + parser=inputParser; + parser.KeepUnmatched=true; + addRequired(parser,'x',@(x) isreal(x)); + addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); + addRequired(parser,'m',@(x) isreal(x) && isrow(x)); + addRequired(parser,'delta',@(x) isreal(x) && isrow(x)); + addRequired(parser,'sigma',@(x) isreal(x) && isscalar(x)); + addRequired(parser,'c',@(x) isreal(x) && isscalar(x)); + addOptional(parser,'side','lower',@(x) strcmpi(x,'lower') || strcmpi(x,'upper') ); + addParameter(parser,'AbsTol',1e-10,@(x) isreal(x) && isscalar(x) && (x>=0)); + addParameter(parser,'RelTol',1e-6,@(x) isreal(x) && isscalar(x) && (x>=0)); + + parse(parser,x,lambda,m,delta,sigma,c,varargin{:}); + side=parser.Results.side; + AbsTol=parser.Results.AbsTol; + RelTol=parser.Results.RelTol; + + u=[]; % pre-allocate in static workspace + + % define the integrand (lambda, m, delta must be column vectors here) + function f=davies_integrand(u,x,lambda,m,delta,sigma) + theta=sum(m.*atan(lambda*u)+(delta.*(lambda*u))./(1+lambda.^2*u.^2),1)/2-u*x/2; + rho=prod(((1+lambda.^2*u.^2).^(m/4)).*exp(((lambda.^2*u.^2).*delta)./(2*(1+lambda.^2*u.^2))),1) .* exp(u.^2*sigma^2/8); + f=sin(theta)./(u.*rho); + end + + % compute the integral + if any(strcmpi(parser.UsingDefaults,'AbsTol')) && any(strcmpi(parser.UsingDefaults,'RelTol')) + davies_integral=arrayfun(@(x) integral(@(u) davies_integrand(u,x-c,lambda',m',delta',sigma),0,inf),x); + if strcmpi(side,'lower') + p=0.5-davies_integral/pi; + elseif strcmpi(side,'upper') + p=0.5+davies_integral/pi; + end + else + syms u + davies_integral=arrayfun(@(x) vpaintegral(@(u) davies_integrand(u,x-c,lambda',m',delta',sigma),... + u,0,inf,'AbsTol',AbsTol,'RelTol',RelTol,'MaxFunctionCalls',inf),x); + if strcmpi(side,'lower') + p=double(0.5-davies_integral/pi); + elseif strcmpi(side,'upper') + p=double(0.5+davies_integral/pi); + end + end + + flag = p<0 | p>1; + p=max(p,0); + p=min(p,1); + +end \ No newline at end of file diff --git a/gx2cdf_imhof.m b/gx2cdf_imhof.m new file mode 100644 index 0000000..031f60f --- /dev/null +++ b/gx2cdf_imhof.m @@ -0,0 +1,125 @@ +function [p,flag]=gx2cdf_imhof(x,lambda,m,delta,c,varargin) + + % GX2CDF_IMHOF Returns the cdf of a generalized chi-squared (a weighted sum of + % non-central chi-squares), using Imhof's [1961] method. + % + % Abhranil Das + % Center for Perceptual Systems, University of Texas at Austin + % If you use this code, please cite: + % A method to integrate and classify normal distributions. + % + % Usage: + % p=gx2cdf_imhof(x,lambda,m,delta,c) + % p=gx2cdf_imhof(x,lambda,m,delta,c,'upper') + % p=gx2cdf_imhof(x,lambda,m,delta,c,'AbsTol',0,'RelTol',1e-7) + % p=gx2cdf_imhof(x,lambda,m,delta,c,'upper','approx','tail') + % + % Example: + % p=gx2cdf_imhof(25,[1 -5 2],[1 2 3],[2 3 7],0) + % + % Required inputs: + % x points at which to evaluate the cdf + % lambda row vector of coefficients of the non-central chi-squares + % m row vector of degrees of freedom of the non-central chi-squares + % delta row vector of non-centrality paramaters (sum of squares of + % means) of the non-central chi-squares + % c constant term + % + % Optional positional input: + % 'upper' more accurate estimate of the complementary CDF when it's small + % + % Optional name-value inputs: + % 'AbsTol' absolute error tolerance for the output + % 'RelTol' relative error tolerance for the output + % The absolute OR the relative tolerance is satisfied. + % 'approx' set to 'tail' for Pearson's approximation of the tail + % probabilities. Works best for the upper (lower) tail when all + % lambda are positive (negative). + % + % Outputs: + % p computed cdf + % flag =true if output was too close to 0 or 1 to compute exactly with + % default settings. Try stricter tolerances or tail approx. for + % more accuracy. + % + % See also: + % Interactive demos + % gx2cdf_davies + % gx2cdf_ruben + % gx2cdf + + parser = inputParser; + addRequired(parser,'x',@(x) isreal(x)); + addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); + addRequired(parser,'m',@(x) isreal(x) && isrow(x)); + addRequired(parser,'delta',@(x) isreal(x) && isrow(x)); + addRequired(parser,'c',@(x) isreal(x) && isscalar(x)); + addOptional(parser,'side','lower',@(x) strcmpi(x,'lower') || strcmpi(x,'upper') ); + addParameter(parser,'AbsTol',1e-10,@(x) isreal(x) && isscalar(x) && (x>=0)); + addParameter(parser,'RelTol',1e-6,@(x) isreal(x) && isscalar(x) && (x>=0)); + addParameter(parser,'approx','none',@(x) strcmpi(x,'none') || strcmpi(x,'tail')); + + parse(parser,x,lambda,m,delta,c,varargin{:}); + side=parser.Results.side; + AbsTol=parser.Results.AbsTol; + RelTol=parser.Results.RelTol; + approx=parser.Results.approx; + + u=[]; % pre-allocate in static workspace + + % define the integrand (lambda, m, delta must be column vectors here) + function f=imhof_integrand(u,x,lambda,m,delta) + theta=sum(m.*atan(lambda*u)+(delta.*(lambda*u))./(1+lambda.^2*u.^2),1)/2-u*x/2; + rho=prod(((1+lambda.^2*u.^2).^(m/4)).*exp(((lambda.^2*u.^2).*delta)./(2*(1+lambda.^2*u.^2))),1); + f=sin(theta)./(u.*rho); + end + + if strcmpi(approx,'tail') % compute tail approximations + j=(1:3)'; + k=sum((lambda.^j).*(j.*delta+m),2); + h=k(2)^3/k(3)^2; + if k(3)>0 + y=(x-c-k(1))*sqrt(h/k(2))+h; + if strcmpi(side,'lower') + p=chi2cdf(y,h); + elseif strcmpi(side,'upper') + p=chi2cdf(y,h,'upper'); + end + else + k=sum(((-lambda).^j).*(j.*delta+m),2); + y=(-(x-c)-k(1))*sqrt(h/k(2))+h; + if strcmpi(side,'lower') + p=chi2cdf(y,h,'upper'); + elseif strcmpi(side,'upper') + p=chi2cdf(y,h); + end + end + + else + % compute the integral + if any(strcmpi(parser.UsingDefaults,'AbsTol')) && any(strcmpi(parser.UsingDefaults,'RelTol')) + imhof_integral=arrayfun(@(x) integral(@(u) imhof_integrand(u,x-c,lambda',m',delta'),0,inf),x); + if strcmpi(side,'lower') + p=0.5-imhof_integral/pi; + elseif strcmpi(side,'upper') + p=0.5+imhof_integral/pi; + end + else + syms u + imhof_integral=arrayfun(@(x) vpaintegral(@(u) imhof_integrand(u,x-c,lambda',m',delta'),... + u,0,inf,'AbsTol',AbsTol,'RelTol',RelTol,'MaxFunctionCalls',inf),x); + + if strcmpi(side,'lower') + p=double(0.5-imhof_integral/pi); + elseif strcmpi(side,'upper') + p=double(0.5+imhof_integral/pi); + end + end + end + + flag = p<0 | p>1; + p=max(p,0); + p=min(p,1); + +end \ No newline at end of file diff --git a/gx2cdf_ruben.m b/gx2cdf_ruben.m new file mode 100644 index 0000000..9a7fbc0 --- /dev/null +++ b/gx2cdf_ruben.m @@ -0,0 +1,90 @@ +function [p,errbnd]=gx2cdf_ruben(x,lambda,m,delta,c,varargin) + + % GX2CDF_RUBEN Returns the cdf of a generalized chi-squared (a weighted sum of + % non-central chi-squares with all weights the same sign), using Ruben's + % [1962] method. + % + % Abhranil Das + % Center for Perceptual Systems, University of Texas at Austin + % If you use this code, please cite: + % A method to integrate and classify normal distributions. + % + % Usage: + % p=gx2cdf_ruben(x,lambda,m,delta,c) + % p=gx2cdf_ruben(x,lambda,m,delta,c,N) + % [p,err]=gx2cdf_ruben(x,lambda,m,delta,c) + % + % Example: + % [p,err]=gx2cdf_ruben(25,[1 5 2],[1 2 3],[2 3 7],0,100) + % + % Required inputs: + % x points at which to evaluate the cdf + % lambda row vector of coefficients of the non-central chi-squares + % m row vector of degrees of freedom of the non-central chi-squares + % delta row vector of non-centrality paramaters (sum of squares of + % means) of the non-central chi-squares + % c constant term + % + % Optional positional input: + % N no. of terms in the approximation. Default = 1000. + % + % Outputs: + % p computed cdf + % errbnd upper error bound of the computed cdf + % + % See also: + % Interactive demos + % gx2cdf_davies + % gx2cdf_imhof + % gx2cdf + + parser=inputParser; + parser.KeepUnmatched=true; + addRequired(parser,'x',@(x) isreal(x)); + addRequired(parser,'lambda',@(x) isreal(x) && isrow(x) && (all(x>0)||all(x<0)) ); + addRequired(parser,'m',@(x) isreal(x) && isrow(x)); + addRequired(parser,'delta',@(x) isreal(x) && isrow(x)); + addRequired(parser,'c',@(x) isreal(x) && isscalar(x)); + addOptional(parser,'side','lower',@(x) strcmpi(x,'lower') || strcmpi(x,'upper') ); + addParameter(parser,'N',1e2,@(x) ismember(x,1:x)); + + parse(parser,x,lambda,m,delta,c,varargin{:}); + side=parser.Results.side; + N=parser.Results.N; + lambda_pos=true; + + if all(lambda<0) + lambda=-lambda; x=-x; lambda_pos=false; + end + beta=0.90625*min(lambda); + M=sum(m); + n=(1:N-1)'; + + % compute the g's + g=sum(m.*(1-beta./lambda).^n,2)+ beta*n.*((1-beta./lambda).^(n-1))*(delta./lambda)'; + + % compute the expansion coefficients + a=nan(N,1); + a(1)=sqrt(exp(-sum(delta))*beta^M*prod(lambda.^(-m))); + if a(1) + % Center for Perceptual Systems, University of Texas at Austin + % If you use this code, please cite: + % A method to integrate and classify normal distributions. + % + % Usage: + % x=gx2inv(p,lambda,m,delta,sigma,c) + % x=gx2inv(p,lambda,m,delta,sigma,c,'AbsTol',0,'RelTol',1e-7) + % + % Example: + % x=gx2inv(0.9,[1 -5 2],[1 2 3],[2 3 7],5,0) + % + % Required inputs: + % p probabilities at which to evaluate the inverse cdf + % lambda row vector of coefficients of the non-central chi-squares + % m row vector of degrees of freedom of the non-central chi-squares + % delta row vector of non-centrality paramaters (sum of squares of + % means) of the non-central chi-squares + % sigma sd of normal term + % c constant term + % + % Optional name-value inputs: + % 'AbsTol' absolute error tolerance for the cdf function that is inverted + % 'RelTol' relative error tolerance for the cdf function that is inverted + % The absolute OR the relative tolerance is satisfied. + % + % Output: + % x computed point + % + % See also: + % Interactive demos + % gx2cdf + % gx2pdf + % gx2rnd + % gx2stat + + parser=inputParser; + parser.KeepUnmatched=true; + addRequired(parser,'p',@(x) isreal(x) && all(x>=0 & x<=1)); + addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); + addRequired(parser,'m',@(x) isreal(x) && isrow(x)); + addRequired(parser,'delta',@(x) isreal(x) && isrow(x)); + addRequired(parser,'sigma',@(x) isreal(x) && isscalar(x)); + addRequired(parser,'c',@(x) isreal(x) && isscalar(x)); + + parse(parser,p,lambda,m,delta,sigma,c,varargin{:}); + + if ~sigma && length(unique(lambda))==1 + % native ncx2 fallback + if sign(unique(lambda))==1 + x=ncx2inv(p,sum(m),sum(delta))*unique(lambda)+c; + elseif sign(unique(lambda))==-1 + x=ncx2inv(1-p,sum(m),sum(delta))*unique(lambda)+c; + end + else + mu=gx2stat(lambda,m,delta,sigma,c); + x=arrayfun(@(p) fzero(@(x) gx2cdf(x,lambda,m,delta,sigma,c,varargin{:})-p,mu),p); + end \ No newline at end of file diff --git a/gx2pdf.m b/gx2pdf.m new file mode 100644 index 0000000..30f48cc --- /dev/null +++ b/gx2pdf.m @@ -0,0 +1,66 @@ +function f=gx2pdf(x,lambda,m,delta,sigma,c,varargin) + + % GX2PDF Returns the pdf of a generalized chi-squared (a weighted sum of + % non-central chi-squares and a normal). + % + % Abhranil Das + % Center for Perceptual Systems, University of Texas at Austin + % If you use this code, please cite: + % A method to integrate and classify normal distributions. + % + % Usage: + % f=gx2pdf(x,lambda,m,delta,sigma,c) + % f=gx2pdf(x,lambda,m,delta,sigma,c,'dx',1e-3) + % f=gx2pdf(x,lambda,m,delta,sigma,c,'AbsTol',0,'RelTol',1e-7) + % + % Example: + % f=gx2pdf(25,[1 -5 2],[1 2 3],[2 3 7],5,0) + % + % Required inputs: + % x points at which to evaluate the pdf + % lambda row vector of coefficients of the non-central chi-squares + % m row vector of degrees of freedom of the non-central chi-squares + % delta row vector of non-centrality paramaters (sum of squares of + % means) of the non-central chi-squares + % sigma sd of normal term + % c constant term + % + % Optional name-value inputs: + % dx step-size for numerically differentiating cdf + % 'AbsTol' absolute error tolerance for the output + % 'RelTol' relative error tolerance for the output + % The absolute OR the relative tolerance is satisfied. + % + % Output: + % f computed pdf + % See also: + % Interactive demos + % gx2cdf + % gx2cdf_davies + % gx2cdf_imhof + % gx2cdf_ruben + + parser = inputParser; + addRequired(parser,'x',@(x) isreal(x)); + addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); + addRequired(parser,'m',@(x) isreal(x) && isrow(x)); + addRequired(parser,'delta',@(x) isreal(x) && isrow(x)); + addRequired(parser,'sigma',@(x) isreal(x) && isscalar(x)); + addRequired(parser,'c',@(x) isreal(x) && isscalar(x)); + addParameter(parser,'AbsTol',1e-10,@(x) isreal(x) && isscalar(x) && (x>=0)); + addParameter(parser,'RelTol',1e-6,@(x) isreal(x) && isscalar(x) && (x>=0)); + [~,v]=gx2stat(lambda,m,delta,sigma,c); + addParameter(parser,'dx',sqrt(v)/100,@(x) isreal(x) && isscalar(x) && (x>=0)); % default derivative step-size is sd/100. + parse(parser,x,lambda,m,delta,sigma,c,varargin{:}); + + if ~sigma && length(unique(lambda))==1 + f=ncx2pdf((x-c)/unique(lambda),sum(m),sum(delta))/abs(unique(lambda)); + else + dx=parser.Results.dx; + p_left=gx2cdf(x-dx,lambda,m,delta,sigma,c,varargin{:}); + p_right=gx2cdf(x+dx,lambda,m,delta,sigma,c,varargin{:}); + f=max((p_right-p_left)/(2*dx),0); + end + +end \ No newline at end of file diff --git a/gx2rnd.m b/gx2rnd.m new file mode 100644 index 0000000..d4895d9 --- /dev/null +++ b/gx2rnd.m @@ -0,0 +1,55 @@ +function r=gx2rnd(lambda,m,delta,sigma,c,varargin) + + % GX2RND Generates generalized chi-squared random numbers. + % + % Abhranil Das + % Center for Perceptual Systems, University of Texas at Austin + % If you use this code, please cite: + % A method to integrate and classify normal distributions. + % + % Usage: + % r=gx2rnd(lambda,m,delta,sigma,c) + % r=gx2rnd(lambda,m,delta,sigma,c,sz) + % r=gx2rnd(lambda,m,delta,sigma,c,sz1,sz2,...) + % r=gx2rnd(lambda,m,delta,sigma,c,[sz1,sz2,...]) + % + % Example: + % r=gx2rnd([1 -5 2],[1 2 3],[2 3 7],5,1,5) + % + % Required inputs: + % lambda row vector of coefficients of the non-central chi-squares + % m row vector of degrees of freedom of the non-central chi-squares + % delta row vector of non-centrality paramaters (sum of squares of + % means) of the non-central chi-squares + % sigma sd of normal term + % c constant term + % + % Optional positional input: + % sz size(s) of the requested array + % + % Output: + % r random number(s) + % + % See also: + % Interactive demos + % gx2stat + % gx2_params_norm_quad + + parser=inputParser; + parser.KeepUnmatched=true; + addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); + addRequired(parser,'m',@(x) isreal(x) && isrow(x)); + addRequired(parser,'delta',@(x) isreal(x) && isrow(x)); + addRequired(parser,'sigma',@(x) isreal(x) && isscalar(x)); + addRequired(parser,'c',@(x) isreal(x) && isscalar(x)); + parse(parser,lambda,m,delta,sigma,c); + + ncxs=arrayfun(@(l,m,d) l*ncx2rnd(m,d,varargin{:}),lambda,m,delta,'un',0); + + r=zeros(size(ncxs{1})); + for i=1:length(ncxs) + r=r+ncxs{i}; + end + r=r+normrnd(c,sigma,varargin{:}); + diff --git a/gx2stat.m b/gx2stat.m new file mode 100644 index 0000000..88f98d1 --- /dev/null +++ b/gx2stat.m @@ -0,0 +1,43 @@ +function [mu,v]=gx2stat(lambda,m,delta,sigma,c) + + % GX2STAT Returns the mean and variance of a generalized chi-squared + % variable (a weighted sum of non-central chi-squares). + % + % Abhranil Das + % Center for Perceptual Systems, University of Texas at Austin + % If you use this code, please cite: + % A method to integrate and classify normal distributions. + % + % Usage: + % [mu,v]=gx2stat(lambda,m,delta,sigma,c) + % + % Example: + % [mu,v]=gx2stat([1 -5 2],[1 2 3],[2 3 7],4,0) + % + % Required inputs: + % lambda row vector of coefficients of the non-central chi-squares + % m row vector of degrees of freedom of the non-central chi-squares + % delta row vector of non-centrality paramaters (sum of squares of + % means) of the non-central chi-squares + % sigma sd of normal term + % c constant term + % + % Outputs: + % mu mean + % v variance + % + % See also: + % Interactive demos + % gx2rnd + % gx2_params_norm_quad + + parser = inputParser; + addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); + addRequired(parser,'m',@(x) isreal(x) && isrow(x)); + addRequired(parser,'delta',@(x) isreal(x) && isrow(x)); + addRequired(parser,'c',@(x) isreal(x) && isscalar(x)); + parse(parser,lambda,m,delta,c); + + mu=dot(lambda,m+delta)+c; + v=2*dot(lambda.^2,m+2*delta)+sigma^2; \ No newline at end of file