diff --git a/Generalized chi-square distribution.mltbx b/Generalized chi-square distribution.mltbx index 422e8f4..841e4db 100644 Binary files a/Generalized chi-square distribution.mltbx and b/Generalized chi-square distribution.mltbx differ diff --git a/Phibar_ray_split.m b/Phibar_ray_split.m new file mode 100644 index 0000000..df7a521 --- /dev/null +++ b/Phibar_ray_split.m @@ -0,0 +1,19 @@ +function [f_big,f_small]=Phibar_ray_split(z,dim) + % returns the complementary cdf of the standard multinormal on a ray, split into a big chunk, which is + % 0, +/-2 or +/-1, and a small part, which is a tail CDF. This is to + % prevent the small part vanishing when the two are added together due to machine imprecision. + + z_c=chi2inv(0.5,dim); % above this criterion, we'll use chi2cdf('upper'). + f_big=(z<=-z_c)+(z=z_c + + f_small=chi2cdf(z.^2,dim); + f_small_upper=chi2cdf(z.^2,dim,'upper'); + f_small(abs(z)>=z_c)=f_small_upper(abs(z)>=z_c); + f_small=f_small.*(1-2*((z<=-z_c)|((z>0)&(zAbhranil 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. + Compute the statistics, pdf, cdf, inverse cdf, random numbers and characteristic function 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.8.8 + 2.0.0 ${PROJECT_ROOT}\Generalized chi-square distribution.mltbx Statistics and Machine Learning Toolbox @@ -123,24 +123,31 @@ LICENSE ${PROJECT_ROOT}\demos.xml ${PROJECT_ROOT}\doc + ${PROJECT_ROOT}\gx2_ellipse.m + ${PROJECT_ROOT}\gx2_ifft.m + ${PROJECT_ROOT}\gx2_imhof.m + ${PROJECT_ROOT}\gx2_imhof_integrand.m + ${PROJECT_ROOT}\gx2_ray_integrand.m + ${PROJECT_ROOT}\gx2_ruben.m ${PROJECT_ROOT}\gx2_to_norm_quad_params.m ${PROJECT_ROOT}\gx2cdf.m - ${PROJECT_ROOT}\gx2cdf_imhof.m - ${PROJECT_ROOT}\gx2cdf_imhof_integrand.m ${PROJECT_ROOT}\gx2cdf_pearson.m ${PROJECT_ROOT}\gx2cdf_ray.m - ${PROJECT_ROOT}\gx2cdf_ruben.m ${PROJECT_ROOT}\gx2char.m ${PROJECT_ROOT}\gx2inv.m + ${PROJECT_ROOT}\gx2pdf.asv ${PROJECT_ROOT}\gx2pdf.m - ${PROJECT_ROOT}\gx2pdf_conv.m - ${PROJECT_ROOT}\gx2pdf_ifft.m - ${PROJECT_ROOT}\gx2pdf_imhof.m - ${PROJECT_ROOT}\gx2pdf_imhof_integrand.m - ${PROJECT_ROOT}\gx2pdf_ray_fast.m + ${PROJECT_ROOT}\gx2pdf_ray.m ${PROJECT_ROOT}\gx2rnd.m ${PROJECT_ROOT}\gx2stat.m + ${PROJECT_ROOT}\int_norm_ray.m + ${PROJECT_ROOT}\log_gx2cdf.m + ${PROJECT_ROOT}\norm_prob_across_angles.m + ${PROJECT_ROOT}\norm_prob_across_rays.m ${PROJECT_ROOT}\norm_quad_to_gx2_params.m + ${PROJECT_ROOT}\phi_ray.m + ${PROJECT_ROOT}\Phibar_ray_split.m + ${PROJECT_ROOT}\standard_quad.m diff --git a/gx2_ellipse.m b/gx2_ellipse.m new file mode 100644 index 0000000..8f24332 --- /dev/null +++ b/gx2_ellipse.m @@ -0,0 +1,104 @@ +function [p,p_rel_err]=gx2_ellipse(x,w,r,lambda,m,varargin) + +parser=inputParser; +parser.KeepUnmatched=true; +addRequired(parser,'x',@(x) isreal(x)); +addRequired(parser,'w',@(x) isreal(x) && isrow(x) && (all(x>0)||all(x<0)) ); +addRequired(parser,'k',@(x) isreal(x) && isrow(x)); +addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); +addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); +addOptional(parser,'side','lower',@(x) strcmpi(x,'lower') || strcmpi(x,'upper') ); +addParameter(parser,'output','cdf',@(x) strcmpi(x,'cdf') || strcmpi(x,'pdf') ); +addParameter(parser,'x_scale','linear',@(x) strcmpi(x,'linear') || strcmpi(x,'log') ); + +parse(parser,x,w,r,lambda,m,varargin{:}); +side=parser.Results.side; +x_scale=parser.Results.x_scale; +w_pos=true; + +% flatten x: +x_flat=x(:); + +if all(w<0) + w=-w; m=-m; w_pos=false; + if strcmpi(x_scale,'linear') + x_flat=-x_flat; + end +end + +% find ellipse parameters +ellipse_center=arrayfun(@(lambda,k) [sqrt(lambda) zeros(1,k-1)],lambda,r,'un',0); +ellipse_center=[ellipse_center{:}]; +ellipse_weights=arrayfun(@(w,k) w*ones(1,k),w,r,'un',0); +ellipse_weights=[ellipse_weights{:}]; + +dim=sum(r); +% factor corr. to the density of the chosen point: +% if norm(ellipse_center) % if non-central +% point_factor=exp(-norm(ellipse_center)^2/2); +% else % if central +% point_factor=exp(-x_flat/(4*min(ellipse_weights))); +% end + +% common factor to cdf and pdf: +% a=point_factor/(2^(dim/2)*gamma(dim/2+1)*sqrt(prod(ellipse_weights))); +a=exp(-norm(ellipse_center)^2/2)/(2^(dim/2)*gamma(dim/2+1)*sqrt(prod(ellipse_weights))); +% ellipse_vol=x_flat.^(dim/2)*pi^(dim/2)/(gamma(dim/2+1)*sqrt(prod(ellipse_weights))); + +if strcmpi(x_scale,'linear') + if strcmpi(parser.Results.output,'cdf') + % p=mvnpdf(ellipse_center)*ellipse_vol; + p=a.*(x_flat-m).^(dim/2); + + % flip if necessary + if (w_pos && strcmpi(side,'upper')) || (~w_pos && strcmpi(side,'lower')) + p=1-p; + end + elseif strcmpi(parser.Results.output,'pdf') + p=(a*dim/2)*(x_flat-m).^(dim/2-1); + end + + % compute log p for log x +else + log10_x=x_flat; + % factor corr. to the density of the chosen point: + % if norm(ellipse_center) % if non-central + % point_factor=norm(ellipse_center)^2; + % else % if central + % point_factor=-10.^log10_x/(2*min(ellipse_weights)); + % end + + if strcmpi(parser.Results.output,'cdf') + log10_p=dim/2*(log10_x-log10(2)) - norm(ellipse_center)^2/log(100) - log10(gamma(dim/2+1)) - sum(log10(ellipse_weights))/2; + elseif strcmpi(parser.Results.output,'pdf') + log10_p=(dim/2-1)*log10_x-(dim/2+1)*log10(2) +log10(dim) - norm(ellipse_center)^2/log(100) - log10(gamma(dim/2+1)) - sum(log10(ellipse_weights))/2; + end + p=log10_p; +end + +% compute the lower and upper bounds for the estimate + +if norm(ellipse_center) % if non-central + if strcmpi(x_scale,'linear') + r=sqrt(x_flat-m)/sqrt(sum(ellipse_center.^2.*ellipse_weights)); + p_rel_err=norm(ellipse_center)^2*r; + % compute log (rel err) for log x + else + p_rel_err=2*log10(norm(ellipse_center))+log10_x/2-log10(sum(ellipse_center.^2.*ellipse_weights))/2; + end +else % if central + p_rel_err=(x_flat-m)/2*min(ellipse_weights); + % compute log (rel err) for log x + if strcmpi(x_scale,'log') + p_rel_err=log10_x-log10(2*min(ellipse_weights)); + end +end + +% p_rel_err(x_flat<0)=log10_p_rel_err; + +% p_rel_err_neg=1-exp((norm(ellipse_center)^2-vecnorm(rect_far,2,2).^2)/2); +% p_rel_err_pos=exp((norm(ellipse_center)^2-vecnorm(rect_near,2,2).^2)/2)-1; + +% reshape outputs to input shape +p=reshape(p,size(x)); +p_rel_err=reshape(p_rel_err,size(x)); \ No newline at end of file diff --git a/gx2_ifft.m b/gx2_ifft.m new file mode 100644 index 0000000..23c0cf8 --- /dev/null +++ b/gx2_ifft.m @@ -0,0 +1,85 @@ +function [p,xgrid]=gx2_ifft(x,w,k,lambda,s,m,varargin) + +parser=inputParser; +parser.KeepUnmatched=true; +addRequired(parser,'x',@(x) isreal(x) || strcmpi(x,'full')); +addRequired(parser,'w',@(x) isreal(x) && isrow(x)); +addRequired(parser,'k',@(x) isreal(x) && isrow(x)); +addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); +addRequired(parser,'s',@(x) isreal(x) && isscalar(x)); +addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); +addOptional(parser,'side','lower',@(x) strcmpi(x,'lower') || strcmpi(x,'upper') ); +addParameter(parser,'output','cdf',@(x) strcmpi(x,'cdf') || strcmpi(x,'pdf') ); +% range of x over which to compute ifft: +if strcmpi(x,'full') + % default span is upto 100 sd from mean + [mu,v]=gx2stat(w,k,lambda,s,m); + span=max(abs(mu+[-1 1]*100*sqrt(v))); +else + % default span is 100* input span + span=100*max(abs(x(:))); +end +addParameter(parser,'span',span,@(x) isscalar(x) && (x>0)); +% number of grid points for ifft +addParameter(parser,'n_grid',1e6+1,@(x) isscalar(x) && (x>0)); +addParameter(parser,'ft_type','cft'); +parse(parser,x,w,k,lambda,s,m,varargin{:}); +span=parser.Results.span; +n_grid=round(parser.Results.n_grid); + +% make n_grid odd +if ~mod(n_grid,2) + n_grid=n_grid+1; +end + +n=(n_grid-1)/2; +idx=-n:n; +dx=span/n; + +xgrid=idx*dx; % grid of x over which this cdf/pdf is computed + +if strcmpi(parser.Results.ft_type,'dft') + % compute non-central and normal pdf's over this span + w_idx=find(w); % non-zero w indices + ncpdfs=nan(nnz(w),length(xgrid)); + for i=1:nnz(w) + if i==1 % add the offset m to the first term + pdf=gx2pdf(xgrid,w(w_idx(i)),k(w_idx(i)),lambda(w_idx(i)),0,m); + else + pdf=gx2pdf(xgrid,w(w_idx(i)),k(w_idx(i)),lambda(w_idx(i)),0,0); + end + pdf(isinf(pdf))=max(pdf(~isinf(pdf))); + ncpdfs(i,:)=pdf; + end + if s + ncpdfs=[ncpdfs;normpdf(xgrid,0,abs(s))]; + end + phi=prod(fft(ifftshift(ncpdfs,2),[],2),1); + p=fftshift(ifft(phi)); + + p=p/(sum(p)*dx); % normalize +elseif strcmpi(parser.Results.ft_type,'cft') + dt=2*pi/(n_grid*dx); + t=idx*dt; % grid of points to evaluate characteristic function + phi=gx2char(-t,w,k,lambda,s,m); + + if strcmpi(parser.Results.output,'pdf') + p=fftshift(ifft(ifftshift(phi)))/dx; + elseif strcmpi(parser.Results.output,'cdf') + phi=phi./(1i*t); + phi(isinf(phi))=0; + p=0.5+fftshift(ifft(ifftshift(phi)))/dx; + end +end + +if strcmpi(parser.Results.output,'cdf') && strcmpi(parser.Results.side,'upper') + p=1-p; +end + + +if ~strcmpi(x,'full') + F=griddedInterpolant(xgrid,p); + p=F(x); +end + +% f=max(f,0); \ No newline at end of file diff --git a/gx2_imhof.m b/gx2_imhof.m new file mode 100644 index 0000000..9030cce --- /dev/null +++ b/gx2_imhof.m @@ -0,0 +1,51 @@ +function [p,errflag]=gx2_imhof(x,w,k,lambda,s,m,varargin) + + parser=inputParser; + parser.KeepUnmatched=true; + addRequired(parser,'x',@(x) isreal(x)); + addRequired(parser,'w',@(x) isreal(x) && isrow(x)); + addRequired(parser,'k',@(x) isreal(x) && isrow(x)); + addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); + addRequired(parser,'s',@(x) isreal(x) && isscalar(x)); + addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); + addOptional(parser,'side','lower',@(x) strcmpi(x,'lower') || strcmpi(x,'upper') ); + addParameter(parser,'output','cdf',@(x) strcmpi(x,'cdf') || strcmpi(x,'pdf') ); + addParameter(parser,'vpa',false,@islogical); + addParameter(parser,'AbsTol',1e-10,@(x) isreal(x) && isscalar(x) && (x>=0)); + addParameter(parser,'RelTol',1e-2,@(x) isreal(x) && isscalar(x) && (x>=0)); + + parse(parser,x,w,k,lambda,s,m,varargin{:}); + output=parser.Results.output; + side=parser.Results.side; + vpaflag=parser.Results.vpa; + AbsTol=parser.Results.AbsTol; + RelTol=parser.Results.RelTol; + + % compute the integral + if ~vpaflag + imhof_integral=arrayfun(@(x) integral(@(u) gx2_imhof_integrand(u,x,w',k',lambda',s,m,output),0,inf,'AbsTol',AbsTol,'RelTol',RelTol),x); + else + syms u + imhof_integral=arrayfun(@(x) vpaintegral(@(u) gx2_imhof_integrand(u,x,w',k',lambda',s,m,output),... + u,0,inf,'AbsTol',AbsTol,'RelTol',RelTol),x); + end + + if strcmpi(output,'cdf') + if strcmpi(side,'lower') + p=double(0.5-imhof_integral/pi); + elseif strcmpi(side,'upper') + p=double(0.5+imhof_integral/pi); + end + errflag = p<0 | p>1; + p=min(p,1); + elseif strcmpi(output,'pdf') + p=imhof_integral/(2*pi); + errflag = p<0; + end + + if any(errflag) + p=max(p,0); + warning('Imhof method output(s) too close to limit to compute exactly, so clipping. Check the flag output, and try stricter tolerances.') + end + +end \ No newline at end of file diff --git a/gx2_imhof_integrand.m b/gx2_imhof_integrand.m new file mode 100644 index 0000000..402fb04 --- /dev/null +++ b/gx2_imhof_integrand.m @@ -0,0 +1,10 @@ +function f=gx2_imhof_integrand(u,x,w,k,lambda,s,m,output) +% define the Imhof integrand (w, k, lambda must be column vectors here) +theta=sum(k.*atan(w*u)+(lambda.*(w*u))./(1+w.^2*u.^2),1)/2+u*(m-x)/2; +rho=prod(((1+w.^2*u.^2).^(k/4)).*exp(((w.^2*u.^2).*lambda)./(2*(1+w.^2*u.^2))),1) .* exp(u.^2*s^2/8); +if strcmpi(output,'cdf') + f=sin(theta)./(u.*rho); +elseif strcmpi(output,'pdf') + f=cos(theta)./rho; +end +end \ No newline at end of file diff --git a/gx2_ray_integrand.m b/gx2_ray_integrand.m new file mode 100644 index 0000000..e315670 --- /dev/null +++ b/gx2_ray_integrand.m @@ -0,0 +1,94 @@ +function [p_rays,p_sym_sum,sym_idx]=gx2_ray_integrand(x,n_z,quad,varargin) +% return the differential probability or probability density on each ray +% that is integrated across rays + +% parse inputs +parser=inputParser; +parser.KeepUnmatched=true; +addRequired(parser,'x',@isnumeric); % points at which to find the cdf/pdf +addRequired(parser,'n_z',@isnumeric); +addRequired(parser,'quad',@isstruct); +addOptional(parser,'side','lower',@(x) strcmpi(x,'lower') || strcmpi(x,'upper') ); +addParameter(parser,'output','prob'); % probability or probability density +addParameter(parser,'vpa',false,@islogical); + +parse(parser,x,n_z,quad,varargin{:}); + +output=parser.Results.output; +side=parser.Results.side; +vpaflag=parser.Results.vpa; + +dim=numel(quad.q1); + +% find the quadratic coefficients across all rays +q2=dot(n_z,quad.q2*n_z,1); +q1=quad.q1'*n_z; +q0=quad.q0; + +% discriminant of the quadratic across rays and levels +delta2=q1.^2-4*q2.*(q0-x); % delta^2 +root_exists=delta2>0; % levels where linear or quadratic roots exist +quad_root_exists=root_exists & q2; % levels where quadratic roots exist +delta=nan(size(delta2)); +delta(quad_root_exists)=sqrt(delta2(quad_root_exists)); % populate only with quad delta for now + +% linear_root_exists=repmat(~q2 & q1, [numel(x) 1]); % levels where linear roots exist + +% sorted roots across rays +z=(-q1+cat(3,-1,1).*delta)./(2*abs(q2)); % quadratic roots where q2 ~= 0 +if nnz(~q2) + z(:,~q2,1)=-(q0-x)./q1(~q2); % linear roots where q2=0 +end + +if strcmpi(output,'prob') + init_sign_rays=sign(4*sign(q2)-2*sign(q1)+sign(q0-x)); + [f_big,f_small]=Phibar_ray_split(z,dim); + p_rays_big=init_sign_rays+1+init_sign_rays.*(f_big(:,:,2)-f_big(:,:,1)); + p_rays_small=init_sign_rays.*(f_small(:,:,2)-f_small(:,:,1)); + if strcmpi(side,'upper') + p_rays=p_rays_big+p_rays_small; + elseif strcmpi(side,'lower') + p_rays=2-p_rays_big-p_rays_small; + end +elseif strcmpi(output,'prob_dens') + + % phi_ray at each root + sum_phi=sum(phi_ray(z,dim),3,'omitnan'); + + % quadratic slope at each root + % quad_slope=nan(numel(x),size(n_z,2)); + % quad_slope(quad_root_exists)=delta(quad_root_exists); + % quad_slope(linear_root_exists)=abs(q1(linear_root_exists)); + quad_slope=nan(size(delta2)); + quad_slope(root_exists)=sqrt(delta2(root_exists)); % slope is the same formula at quad and linear roots + + % divide phi_ray by quad slope + p_rays=sum_phi./quad_slope; + + p_rays(isnan(p_rays))=0; +end + +p_sym_sum=num2cell(zeros(numel(x),1)); +sym_idx=[]; +% cases where root exists but computed prob. is 0, +tiny_probs=root_exists&(~p_rays); +if nnz(tiny_probs) + if ~vpaflag + warning("%.1f%% of rays contain probabilities smaller than realmin=1e-308, returning 0. Set 'vpa' to true to compute these with variable precision.",100*mean(tiny_probs,'all')) + else + sym_idx=any(tiny_probs,2); % levels at which at least one ray is sym. + % group roots into a cell array so we can use cellfun: + z_cell=permute(num2cell(permute(z,[1 3 2]),2),[1 3 2]); + % remove nan roots: + z_cell=cellfun(@(z_each) z_each(~isnan(z_each)),z_cell,'un',0); + if strcmpi(output,'prob') + p_sym_sum(sym_idx)=sym2cell(sum(cellfun(@(init_sign_ray,z_ray) prob_ray_sym(init_sign_ray,z_ray,dim,side),num2cell(init_sign_rays(sym_idx,:)),z_cell(sym_idx,:)),2)); + elseif strcmpi(output,'prob_dens') + % p_sym_sum(sym_idx)=sym2cell(arrayfun(@(iLevel) ... + % sum(sum(phi_ray(sym(z(iLevel,tiny_probs(iLevel,:),:)),dim),3)./quad_slope(iLevel,tiny_probs(iLevel,:))),... + % find(sym_idx))); + p_sym=cellfun(@(z_ray,quad_slope_ray) sum(phi_ray(sym(z_ray),dim))/quad_slope_ray,z_cell(sym_idx,:),num2cell(quad_slope(sym_idx,:))); + p_sym_sum(sym_idx)=sym2cell(sum(p_sym,2,'omitnan')); + end + end +end \ No newline at end of file diff --git a/gx2_ruben.m b/gx2_ruben.m new file mode 100644 index 0000000..0b3e98f --- /dev/null +++ b/gx2_ruben.m @@ -0,0 +1,75 @@ +function [p,p_err]=gx2_ruben(x,w,k,lambda,m,varargin) + +parser=inputParser; +parser.KeepUnmatched=true; +addRequired(parser,'x',@(x) isreal(x)); +addRequired(parser,'w',@(x) isreal(x) && isrow(x) && (all(x>0)||all(x<0)) ); +addRequired(parser,'k',@(x) isreal(x) && isrow(x)); +addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); +addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); +addOptional(parser,'side','lower',@(x) strcmpi(x,'lower') || strcmpi(x,'upper') ); +addParameter(parser,'output','cdf',@(x) strcmpi(x,'cdf') || strcmpi(x,'pdf') ); +addParameter(parser,'n_ruben',1e2,@(x) ismember(x,1:x)); + +parse(parser,x,w,k,lambda,m,varargin{:}); +side=parser.Results.side; +n_ruben=parser.Results.n_ruben; + +% flatten x: +x_flat=x(:); + +w_pos=true; +if all(w<0) + w=-w; x_flat=-x_flat; m=-m; w_pos=false; +end +beta=0.90625*min(w); +M=sum(k); +n=(1:n_ruben-1)'; + +% compute the g's +g=sum(k.*(1-beta./w).^n,2)+ beta*n.*((1-beta./w).^(n-1))*(lambda./w)'; + +% compute the expansion coefficients +a=nan(n_ruben,1); +a(1)=sqrt(exp(-sum(lambda))*beta^M*prod(w.^(-k))); +if a(1) + % Abhranil Das % Center for Perceptual Systems, University of Texas at Austin + % Comments, questions, bugs to abhranil.das@utexas.edu % If you use this code, please cite: - % A method to integrate and classify normal distributions. + % 1. A method to integrate and classify normal distributions + % 2. New methods for computing the generalized chi-square distribution % - % Usage: - % quad=gx2_to_norm_quad_params(w,k,lambda,m,s) + % Usage: + % quad=gx2_to_norm_quad_params(w,k,lambda,s,m) % % Example: - % quad=gx2_to_norm_quad_params([1 -5 2],[1 2 3],[2 3 7],5,10); + % quad=gx2_to_norm_quad_params([1 -5 2],[1 2 3],[2 3 7],10,5); + % + % Required inputs: + % w row vector of weights of the non-central chi-squares + % k row vector of degrees of freedom of the non-central chi-squares + % lambda row vector of non-centrality paramaters (sum of squares of + % means) of the non-central chi-squares + % s scale of normal term + % m offset % - % Required inputs: - % w row vector of weights of the non-central chi-squares - % k row vector of degrees of freedom of the non-central chi-squares - % lambda row vector of non-centrality paramaters (sum of squares of - % means) of the non-central chi-squares - % m mean of normal term - % s sd of normal term - % % % Outputs: - % quad struct with quadratic form coefficients: - % q2 matrix of quadratic coefficients - % q1 column vector of linear coefficients - % q0 scalar constant + % quad struct with quadratic form coefficients: + % q2 matrix of quadratic coefficients + % q1 column vector of linear coefficients + % q0 scalar constant % The dimension of the standard normal is the length of q1. % % See also: - % Interactive demos - % gx2rnd - % gx2stat - % gx2cdf - % gx2pdf + % Getting Started guide q2_parts=arrayfun(@(w,k) w*ones(1,k),w,k,'un',0); % each w_i, k_i times q2=horzcat(q2_parts{:}); % append all of them diff --git a/gx2cdf.m b/gx2cdf.m index bc7fcc7..dd35383 100644 --- a/gx2cdf.m +++ b/gx2cdf.m @@ -1,88 +1,147 @@ -function p=gx2cdf(x,w,k,lambda,m,s,varargin) +function [p,p_err,x_grid]=gx2cdf(x,w,k,lambda,s,m,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. + % GX2CDF Returns the cdf of a generalized chi-squared distribution. % - % Abhranil Das + % Abhranil Das % Center for Perceptual Systems, University of Texas at Austin + % Comments, questions, bugs to abhranil.das@utexas.edu % If you use this code, please cite: - % A method to integrate and classify normal distributions. + % 1. A method to integrate and classify normal distributions + % 2. New methods for computing the generalized chi-square distribution % % Usage: - % p=gx2cdf(x,w,k,lambda,m,s) - % p=gx2cdf(x,w,k,lambda,m,s,'upper') - % p=gx2cdf(x,w,k,lambda,m,s,'AbsTol',0,'RelTol',1e-7) + % p=gx2cdf(x,w,k,lambda,s,m) + % p=gx2cdf(x,w,k,lambda,s,m,'upper') + % p=gx2cdf(x,w,k,lambda,s,m,'method','imhof','AbsTol',0,'RelTol',1e-7) + % [p,p_err]=gx2cdf(x,w,k,lambda,s,m,'method','ray','n_rays',1e7) + % [p,~,x_grid]=gx2cdf('full',w,k,lambda,s,m) + % etc. % % Example: - % f=gx2cdf(25,[1 -5 2],[1 2 3],[2 3 7],5,0) + % p=gx2cdf(25,[1 -5 2],[1 2 3],[2 3 7],0,5) % % Required inputs: - % x points at which to evaluate the CDF + % x array of points at which to evaluate the CDF. + % 'full' to use IFFT method to return p over an array of x that spans the distribution. + % if method is 'ellipse' and 'x_scale' is 'log', these are log10 + % values of points measured from the finite tail, i.e. + % log10(abs(x-m)), and will return log10 of p. + % % w row vector of weights of the non-central chi-squares % k row vector of degrees of freedom of the non-central chi-squares % lambda row vector of non-centrality paramaters (sum of squares of % means) of the non-central chi-squares - % m mean of normal term - % s sd of normal term + % s scale of normal term + % m offset % % 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 + % method 'auto' (default) tries to pick the best method for the parameters + % 'imhof' for Imhof-Davies method, works for all parameters + % 'ray' for ray-trace method, works for all parameters + % 'ifft' for IFFT method, works for all parameters + % 'ruben' for Ruben's method. All w must be same sign and s=0. + % 'ellipse' for ellipse approximation. All w must be same sign and s=0. + % vpa true to use variable precision in Imhof and ray methods. Default=false. + % AbsTol absolute error tolerance for the output. Default=1e-10. + % RelTol relative error tolerance for the output. Default=1e-6. + % AbsTol and RelTol are used only for Imhof's method, and ray method with grid integration. % The absolute OR the relative tolerance is satisfied. - % vpa false (default) to do the integral in Imhof's method numerically, - % true to do it symbolically with variable precision. % - % Output: - % p computed cdf + % Options for ray method: + % force_mc true to force Monte-Carlo instead of grid integration over + % rays. Default=false. + % n_rays no. of rays for Monte-Carlo integration. Larger=more accurate. + % gpu_batch no. of rays to compute on each GPU batch. 0 to use CPU. + % + % Options for Ruben method: + % n_ruben no. of terms in Ruben's series. Larger=more accurate. Default=100. + % + % Options for ifft method: + % span span of x over which to compute the IFFT. Larger=more accurate. + % n_grid number of grid points used for IFFT. Larger=more accurate. Default=1e6. + % + % Options for ellipse method: + % x_scale 'linear' (default). 'log' if input x is log10 values of x, to compute on small x values. + % + % Outputs: + % p computed cdf. + % (if ray method with vpa) can be symbolic for small values + % (if ellipse method with 'log' x_scale) log10 of p + % p_err (if ray method with Monte-Carlo integration) standard error of the output p + % (if Ruben's method) upper error bound of the output p + % (if Imhof's method) logical array, true for outputs too close to 0 or 1 to compute exactly with + % default settings. Try stricter tolerances. + % (if ellipse method) relative error of p + % (if ellipse method with 'log' x_scale) log10 of relative error of p + % x_grid if input x is 'full', this returns the x points at which p is returned % % See also: - % Interactive demos - % gx2cdf_davies - % gx2cdf_imhof - % gx2cdf_ruben - % gx2pdf + % Getting Started guide parser = inputParser; parser.KeepUnmatched = true; - addRequired(parser,'x',@(x) isreal(x)); - addRequired(parser,'w',@(x) isreal(x) && isrow(x)); - addRequired(parser,'k',@(x) isreal(x) && isrow(x)); - addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); - addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); + addRequired(parser,'x',@(x) isreal(x) || strcmpi(x,'full')); + addRequired(parser,'w',@(x) isreal(x)); + addRequired(parser,'k',@(x) isreal(x)); + addRequired(parser,'lambda',@(x) isreal(x)); addRequired(parser,'s',@(x) isreal(x) && isscalar(x)); + addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); addOptional(parser,'side','lower',@(x) strcmpi(x,'lower') || strcmpi(x,'upper') ); addParameter(parser,'method','auto'); - 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,w,k,lambda,m,s,varargin{:}); + parse(parser,x,w,k,lambda,s,m,varargin{:}); method=parser.Results.method; side=parser.Results.side; + if strcmpi(x,'full') + method='ifft'; + end + if strcmpi(method,'auto') - if ~s && length(unique(w))==1 + if ~s && length(unique(w))==1 % no s and only one unique weight % ncx2 fallback if (sign(unique(w))==1 && strcmpi(side,'lower')) || (sign(unique(w))==-1 && strcmpi(side,'upper')) p=ncx2cdf((x-m)/unique(w),sum(k),sum(lambda)); else p=ncx2cdf((x-m)/unique(w),sum(k),sum(lambda),'upper'); end + elseif sum(w)==0 && s % only normal term + if strcmpi(side,'lower') + p=normcdf(x,m,s); + elseif strcmpi(side,'upper') + p=normcdf(x,m,s,'upper'); + end + elseif ~s + if (all(w>0) && strcmpi(side,'lower'))||(all(w<0) && strcmpi(side,'upper')) % no s and w same sign + [p,p_err]=gx2_ruben(x,w,k,lambda,m,varargin{:}); + else + [p,p_err]=gx2_imhof(x,w,k,lambda,s,m,varargin{:}); + end else - p=gx2cdf_ray(x,w,k,lambda,m,s,varargin{:}); + [p,p_err]=gx2_imhof(x,w,k,lambda,s,m,varargin{:}); end + elseif strcmpi(method,'ifft') + [p,x_grid]=gx2_ifft(x,w,k,lambda,s,m,varargin{:},'output','cdf'); + p_err=[]; elseif strcmpi(method,'ray') - p=gx2cdf_ray(x,w,k,lambda,m,s,varargin{:}); + [p,p_err]=gx2cdf_ray(x,w,k,lambda,s,m,varargin{:}); elseif strcmpi(method,'imhof') - p=gx2cdf_imhof(x,w,k,lambda,m,s,varargin{:}); + [p,p_err]=gx2_imhof(x,w,k,lambda,s,m,varargin{:}); elseif strcmpi(method,'ruben') if s || ~(all(w>0)||all(w<0)) error("Ruben's method can only be used when all w are the same sign and s=0.") else - p=gx2cdf_ruben(x,w,k,lambda,m,varargin{:}); + [p,p_err]=gx2_ruben(x,w,k,lambda,m,varargin{:}); + end + elseif strcmpi(method,'ellipse') + if s || ~(all(w>0)||all(w<0)) + error("The ellipse approximation can only be used when all w are the same sign and s=0.") + else + [p,p_err]=gx2_ellipse(x,w,k,lambda,m,varargin{:}); end end \ No newline at end of file diff --git a/gx2cdf_imhof.m b/gx2cdf_imhof.m deleted file mode 100644 index 873e25e..0000000 --- a/gx2cdf_imhof.m +++ /dev/null @@ -1,93 +0,0 @@ -function [p,errflag]=gx2cdf_imhof(x,w,k,lambda,m,s,varargin) - - % GX2CDF_IMHOF Returns the cdf of a generalized chi-squared (a weighted - % sum of non-central chi-squares and a normal), using Davies' [1973] - % method, which extends Imhof's [1961] method to include the s term. - % - % 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,w,k,lambda,m,s) - % p=gx2cdf_imhof(x,w,k,lambda,m,s,'upper') - % p=gx2cdf_imhof(x,w,k,lambda,m,s,'AbsTol',0,'RelTol',1e-7) - % - % Example: - % p=gx2cdf_imhof(25,[1 -5 2],[1 2 3],[2 3 7],5,0) - % - % Required inputs: - % x points at which to evaluate the cdf - % w row vector of weights of the non-central chi-squares - % k row vector of degrees of freedom of the non-central chi-squares - % lambda row vector of non-centrality paramaters (sum of squares of - % means) of the non-central chi-squares - % m mean of normal term - % s sd of normal 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. - % vpa false (default) to do the integral numerically, - % true to do it symbolically with variable precision. - % - % Outputs: - % p computed cdf - % flag = true for output(s) 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,'w',@(x) isreal(x) && isrow(x)); - addRequired(parser,'k',@(x) isreal(x) && isrow(x)); - addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); - addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); - addRequired(parser,'s',@(x) isreal(x) && isscalar(x)); - addOptional(parser,'side','lower',@(x) strcmpi(x,'lower') || strcmpi(x,'upper') ); - addParameter(parser,'vpa',false,@islogical); - 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,w,k,lambda,m,s,varargin{:}); - side=parser.Results.side; - vpaflag=parser.Results.vpa; - AbsTol=parser.Results.AbsTol; - RelTol=parser.Results.RelTol; - - % compute the integral - if ~vpaflag - imhof_integral=arrayfun(@(x) integral(@(u) gx2cdf_imhof_integrand(u,x,w',k',lambda',m,s),0,inf,'AbsTol',AbsTol,'RelTol',RelTol),x); - else - syms u - imhof_integral=arrayfun(@(x) vpaintegral(@(u) gx2cdf_imhof_integrand(u,x,w',k',lambda',m,s),... - u,0,inf,'AbsTol',AbsTol,'RelTol',RelTol,'MaxFunctionCalls',inf),x); - end - - if strcmpi(side,'lower') - p=double(0.5-imhof_integral/pi); - elseif strcmpi(side,'upper') - p=double(0.5+imhof_integral/pi); - end - - errflag = p<0 | p>1; - if any(errflag) - warning('Imhof method gx2cdf output(s) too close to 0/1 to compute exactly, so clipping to 0/1. Check the flag output, and try stricter tolerances.') - p=max(p,0); - p=min(p,1); - end - - -end \ No newline at end of file diff --git a/gx2cdf_imhof_integrand.m b/gx2cdf_imhof_integrand.m deleted file mode 100644 index 0a541e6..0000000 --- a/gx2cdf_imhof_integrand.m +++ /dev/null @@ -1,6 +0,0 @@ -function f=gx2cdf_imhof_integrand(u,x,w,k,lambda,m,s) - % define the Imhof integrand (w, k, lambda must be column vectors here) - theta=sum(k.*atan(w*u)+(lambda.*(w*u))./(1+w.^2*u.^2),1)/2+u*(m-x)/2; - rho=prod(((1+w.^2*u.^2).^(k/4)).*exp(((w.^2*u.^2).*lambda)./(2*(1+w.^2*u.^2))),1) .* exp(u.^2*s^2/8); - f=sin(theta)./(u.*rho); -end \ No newline at end of file diff --git a/gx2cdf_ray.m b/gx2cdf_ray.m index a567ad3..60e5c96 100644 --- a/gx2cdf_ray.m +++ b/gx2cdf_ray.m @@ -1,97 +1,32 @@ -function p=gx2cdf_ray(x,w,k,lambda,m,s,varargin) - - % GX2CDF_RAY Returns the cdf of a generalized chi-squared (a weighted - % sum of non-central chi-squares and a normal), using the ray method [CITE]. - % - % 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,w,k,lambda,m,s) - % p=gx2cdf_davies(x,w,k,lambda,m,s,'upper') - % p=gx2cdf_davies(x,w,k,lambda,m,s,'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 - % w row vector of weights of the non-central chi-squares - % k row vector of degrees of freedom of the non-central chi-squares - % lambda row vector of non-centrality paramaters (sum of squares of - % means) of the non-central chi-squares - % m mean of normal term - % s sd of normal 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 +function [p,p_err]=gx2cdf_ray(x,w,k,lambda,s,m,varargin) parser=inputParser; parser.KeepUnmatched=true; addRequired(parser,'x',@(x) isreal(x)); - addRequired(parser,'w',@(x) isreal(x) && isrow(x)); - addRequired(parser,'k',@(x) isreal(x) && isrow(x)); - addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); - addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); + addRequired(parser,'w',@(x) isreal(x)); + addRequired(parser,'k',@(x) isreal(x)); + addRequired(parser,'lambda',@(x) isreal(x)); addRequired(parser,'s',@(x) isreal(x) && isscalar(x)); + addRequired(parser,'m',@(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-2,@(x) isreal(x) && isscalar(x) && (x>=0)); - addParameter(parser,'mc_samples',500); - - parse(parser,x,w,k,lambda,m,s,varargin{:}); - side=parser.Results.side; - AbsTol=parser.Results.AbsTol; - RelTol=parser.Results.RelTol; - mc_samples=parser.Results.mc_samples; - % find the quadratic form of the standard multinormal - quad=gx2_to_norm_quad_params(w,k,lambda,m,s); - dim=numel(quad.q1); + parse(parser,x,w,k,lambda,s,m,varargin{:}); - % flip the domain, for lower side - if strcmpi(side,'lower') - quad=structfun(@uminus,quad,'un',0); - x=-x; - end + y=x(:); % flatten input array x into a column vector of levels y - % function to compute cdf for each point in x - function p_each=gx2cdf_ray_each(x_each,quad,dim,AbsTol,RelTol,mc_samples) - tic - quad.q0=quad.q0-x_each; - p_each=int_norm_ray(zeros(dim,1),eye(dim),quad,'AbsTol',AbsTol,'RelTol',RelTol,'mc_samples',mc_samples); -% p_each=integrate_normal(zeros(dim,1),eye(dim),quad,'method','ray','AbsTol',AbsTol,'RelTol',RelTol,'mc_samples',mc_samples,'plotmode',0); - toc - end + % find standard quadratic form corresponding to the gx2: + quad=gx2_to_norm_quad_params(w,k,lambda,s,m); + dim=numel(quad.q1); + mu=zeros(dim,1); + v=eye(dim); - % integrate the standard multinormal over the domain using the ray method - p=arrayfun(@(x_each) gx2cdf_ray_each(x_each,quad,dim,AbsTol,RelTol,mc_samples),x,'un',0); + % integrate with 'complement' for lower side. + [p,p_err]=int_norm_ray(mu,v,quad,varargin{:},'fun_level',y); - % if p is a cell but there are no symbols, convert to numeric array - if iscell(p) - num_idx=cellfun(@isnumeric, p); - if all(num_idx) - p=cell2mat(p); - end + % reshape flattened output arrays to shape of input x + p=reshape(p,size(x)); + if ~isempty(p_err) + p_err=reshape(p_err,size(x)); end end \ No newline at end of file diff --git a/gx2cdf_ruben.m b/gx2cdf_ruben.m deleted file mode 100644 index a87c9ab..0000000 --- a/gx2cdf_ruben.m +++ /dev/null @@ -1,90 +0,0 @@ -function [p,errbnd]=gx2cdf_ruben(x,w,k,lambda,m,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,w,k,lambda,m) - % p=gx2cdf_ruben(x,w,k,lambda,m,N) - % [p,err]=gx2cdf_ruben(x,w,k,lambda,m) - % - % 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 - % w row vector of weights of the non-central chi-squares - % k row vector of degrees of freedom of the non-central chi-squares - % lambda row vector of non-centrality paramaters (sum of squares of - % means) of the non-central chi-squares - % m mean of normal 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,'w',@(x) isreal(x) && isrow(x) && (all(x>0)||all(x<0)) ); - addRequired(parser,'k',@(x) isreal(x) && isrow(x)); - addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); - addRequired(parser,'m',@(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,w,k,lambda,m,varargin{:}); - side=parser.Results.side; - N=parser.Results.N; - lambda_pos=true; - - if all(w<0) - w=-w; x=-x; m=-m; lambda_pos=false; - end - beta=0.90625*min(w); - M=sum(k); - n=(1:N-1)'; - - % compute the g's - g=sum(k.*(1-beta./w).^n,2)+ beta*n.*((1-beta./w).^(n-1))*(lambda./w)'; - - % compute the expansion coefficients - a=nan(N,1); - a(1)=sqrt(exp(-sum(lambda))*beta^M*prod(w.^(-k))); - if a(1)A method to integrate and classify normal distributions + % 2. New methods for computing the generalized chi-square distribution + % + % Usage: + % phi=gx2char(t,w,k,lambda,s,m) + % + % Example: + % phi=gx2char(linspace(-10,10,100),[1 -5 2],[1 2 3],[2 3 7],0,5) + % + % Required inputs: + % t array of points at which to compute the characteristic function + % w row vector of weights of the non-central chi-squares + % k row vector of degrees of freedom of the non-central chi-squares + % lambda row vector of non-centrality paramaters (sum of squares of + % means) of the non-central chi-squares + % s scale of normal term + % m offset + % + % Outputs: + % phi characteristic function + % + % See also: + % Getting Started guide parser = inputParser; addRequired(parser,'w',@(x) isreal(x) && isrow(x)); addRequired(parser,'k',@(x) isreal(x) && isrow(x)); addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); - addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); addRequired(parser,'s',@(x) isreal(x) && isscalar(x)); - parse(parser,w,k,lambda,m,s); + addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); + parse(parser,w,k,lambda,s,m); - t_flat=t(:); % flatten input array x into a column vector of levels y + t_flat=t(:); % flatten input array phi=exp(1i*m*t_flat+1i*t_flat.*sum((w.*lambda)./(1-2i*t_flat.*w),2)-s^2*t_flat.^2/2)./... prod((1-2i*w.*t_flat).^(k/2),2); diff --git a/gx2inv.m b/gx2inv.m index 40c2cef..3f9b9a9 100644 --- a/gx2inv.m +++ b/gx2inv.m @@ -1,67 +1,84 @@ -function x=gx2inv(p,w,k,lambda,m,s,varargin) - - % GX2INV Returns the inverse cdf of a generalized chi-squared, using - % Ruben's [1962] method, Davies' [1973] method, or the native ncx2inv, - % 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: - % x=gx2inv(p,w,k,lambda,m,s) - % x=gx2inv(p,w,k,lambda,m,s,'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 - % w row vector of weights of the non-central chi-squares - % k row vector of degrees of freedom of the non-central chi-squares - % lambda row vector of non-centrality paramaters (sum of squares of - % means) of the non-central chi-squares - % m mean of normal term - % s sd of normal 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 +function x=gx2inv(p,w,k,lambda,s,m,varargin) + + % GX2INV Returns the inverse cdf of a generalized chi-squared distribution. + % + % Abhranil Das + % Center for Perceptual Systems, University of Texas at Austin + % Comments, questions, bugs to abhranil.das@utexas.edu + % If you use this code, please cite: + % 1. A method to integrate and classify normal distributions + % 2. New methods for computing the generalized chi-square distribution + % + % Usage: + % x=gx2inv(p,w,k,lambda,s,m) + % x=gx2inv(p,w,k,lambda,s,m,'upper','method','imhof','AbsTol',0,'RelTol',1e-7) + % etc. + % + % Example: + % x=gx2inv(0.9,[1 -5 2],[1 2 3],[2 3 7],0,5) + % x=gx2inv(-100,[1 -5 2],[1 2 3],[2 3 7],0,5,'upper','method','ray','n_rays',1e4) + % + % Required inputs: + % p probabilities at which to evaluate the inverse cdf. + % Negative values indicate log probability. + % w row vector of weights of the non-central chi-squares + % k row vector of degrees of freedom of the non-central chi-squares + % lambda row vector of non-centrality paramaters (sum of squares of + % means) of the non-central chi-squares + % s scale of normal term + % m offset + % + % Optional positional input: + % 'upper' for more accurate quantiles when entering an upper tail + % probability (complementary cdf) + % + % Optional name-value inputs: + % This function numerically finds roots of the gx2cdf function, so most + % options for the gx2cdf function can be used here, eg 'method', which + % will be passed on to gx2cdf + % + % Output: + % x computed quantile % % 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,'w',@(x) isreal(x) && isrow(x)); - addRequired(parser,'k',@(x) isreal(x) && isrow(x)); - addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); - addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); - addRequired(parser,'s',@(x) isreal(x) && isscalar(x)); - - parse(parser,p,w,k,lambda,m,s,varargin{:}); - - if ~s && length(unique(w))==1 - % native ncx2 fallback - if sign(unique(w))==1 - x=ncx2inv(p,sum(k),sum(lambda))*unique(w)+m; - elseif sign(unique(w))==-1 - x=ncx2inv(1-p,sum(k),sum(lambda))*unique(w)+m; + % Getting Started guide + + parser=inputParser; + parser.KeepUnmatched=true; + addRequired(parser,'p',@(x) isreal(x) && all(x<=1)); + addRequired(parser,'w',@(x) isreal(x)); + addRequired(parser,'k',@(x) isreal(x)); + addRequired(parser,'lambda',@(x) isreal(x)); + addRequired(parser,'s',@(x) isreal(x) && isscalar(x)); + addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); + addOptional(parser,'side','lower',@(x) strcmpi(x,'lower') || strcmpi(x,'upper') ); + + parse(parser,p,w,k,lambda,s,m,varargin{:}); + + side=parser.Results.side; + + if ~s && length(unique(w))==1 + % native ncx2 fallback + if strcmpi(side,'upper') + p=1-p; + end + if sign(unique(w))==1 + x=ncx2inv(p,sum(k),sum(lambda))*unique(w)+m; + elseif sign(unique(w))==-1 + x=ncx2inv(1-p,sum(k),sum(lambda))*unique(w)+m; elseif unique(w)==0 x=0; - end - else - mu=gx2stat(w,k,lambda,s,m); - x=arrayfun(@(p) fzero(@(x) gx2cdf(x,w,k,lambda,m,s,varargin{:})-p,mu),p); - end \ No newline at end of file + end + else + mu=gx2stat(w,k,lambda,s,m); + if p>0 + x=arrayfun(@(p) fzero(@(x) gx2cdf(x,w,k,lambda,s,m,varargin{:})-p,mu),p); + else % log probability, inverted using sym + % options = optimset('PlotFcns',{@optimplotx,@optimplotfval}); + % x=arrayfun(@(p) fzero(@(x) double(log10(gx2cdf(x,w,k,lambda,s,m,varargin{:})))-p,mu),p); + % options = optimset('FunValCheck','off') + x=arrayfun(@(p) fzero(@(x) log_gx2cdf(x,w,k,lambda,s,m,varargin{:})-p,mu),p); + end + end \ No newline at end of file diff --git a/gx2pdf.m b/gx2pdf.m index b1e4494..485c778 100644 --- a/gx2pdf.m +++ b/gx2pdf.m @@ -1,87 +1,147 @@ -function [f,xgrid]=gx2pdf(x,w,k,lambda,m,s,varargin) +function [f,f_err,xgrid]=gx2pdf(x,w,k,lambda,s,m,varargin) - % GX2PDF Returns the pdf of a generalized chi-squared (a weighted sum of - % non-central chi-squares and a normal). + % GX2PDF Returns the pdf of a generalized chi-squared distribution. % - % Abhranil Das + % Abhranil Das % Center for Perceptual Systems, University of Texas at Austin + % Comments, questions, bugs to abhranil.das@utexas.edu % If you use this code, please cite: - % A method to integrate and classify normal distributions. + % 1. A method to integrate and classify normal distributions + % 2. New methods for computing the generalized chi-square distribution % % Usage: - % f=gx2pdf(x,w,k,lambda,m,s) - % f=gx2pdf(x,w,k,lambda,m,s,'dx',1e-3) - % [f,xfull]=gx2pdf('full',w,k,lambda,m,s) - % f=gx2pdf(x,w,k,lambda,m,s,'method','conv','AbsTol',0,'RelTol',1e-7) + % f=gx2pdf(x,w,k,lambda,s,m) + % f=gx2pdf(x,w,k,lambda,s,m,'method','imhof','AbsTol',0,'RelTol',1e-7) + % [f,f_err]=gx2pdf(x,w,k,lambda,s,m,'method','ray','n_rays',1e7) + % [f,~,x_grid]=gx2pdf('full',w,k,lambda,s,m) + % etc. % % Example: - % f=gx2pdf(25,[1 -5 2],[1 2 3],[2 3 7],5,0) - % f=gx2pdf([17 25],[1 -5 2],[1 2 3],[2 3 7],5,0,'method','conv') - % [f,xfull]=gx2pdf('full',[1 -5 2],[1 2 3],[2 3 7],5,0) + % f=gx2pdf(25,[1 -5 2],[1 2 3],[2 3 7],0,5) % % Required inputs: - % x points at which to evaluate the pdf, or 'full', to return - % points xfull and f over a full range of the pdf (this uses - % the 'conv' method) + % x array of points at which to evaluate the PDF. + % 'full' to use IFFT method to return f over an array of x that spans the distribution. + % if method is 'ellipse' and 'x_scale' is 'log', these are log10 + % values of points measured from the finite tail, i.e. + % log10(abs(x-m)), and will return log10 of f. + % % w row vector of weights of the non-central chi-squares % k row vector of degrees of freedom of the non-central chi-squares % lambda row vector of non-centrality paramaters (sum of squares of % means) of the non-central chi-squares - % m mean of normal term - % s sd of normal term + % s scale of normal term + % m offset % % Optional name-value inputs: - % method 'diff' (default) for differentiating the generalized - % chi-square cdf, or 'conv' for convolving noncentral - % chi-square pdf's - % dx step-size of fineness (for convolving or differentiating) - % AbsTol absolute error tolerance for the output when using 'diff' - % RelTol relative error tolerance for the output when using 'diff' + % method 'auto' (default) tries to pick the best method for the parameters + % 'imhof' for Imhof-Davies method, works for all parameters + % 'ray' for ray-trace method, works for all parameters + % 'ifft' for IFFT method, works for all parameters + % 'ruben' for Ruben's method. All w must be same sign and s=0. + % 'ellipse' for ellipse approximation. All w must be same sign and s=0. + % diff default=false. True for numerically differentiating the output of gx2cdf, + % in which case options for gx2cdf can be used here, and will + % be passed to gx2cdf. + % dx step-size for numerically differentiating the cdf. Default + % = sd of the distribution divided by 1e4. + % vpa true to use variable precision in Imhof and ray methods. Default=false. + % AbsTol absolute error tolerance for the output. Default=1e-10. + % RelTol relative error tolerance for the output. Default=1e-6. + % AbsTol and RelTol are used only for Imhof's method, and ray method with grid integration. % The absolute OR the relative tolerance is satisfied. % - % Output: - % f computed pdf at points, or over full range - % xfull x-values over the full pdf range, returned when x is 'full' + % Options for ray method: + % force_mc true to force Monte-Carlo instead of grid integration over + % rays. Default=false. + % n_rays no. of rays for Monte-Carlo integration. Larger=more accurate. + % gpu_batch no. of rays to compute on each GPU batch. 0 to use CPU. + % + % Options for Ruben method: + % n_ruben no. of terms in Ruben's series. Larger=more accurate. Default=100. + % + % Options for ifft method: + % span span of x over which to compute the IFFT. Larger=more accurate. + % n_grid number of grid points used for IFFT. Larger=more accurate. Default=1e6. + % + % Options for ellipse method: + % x_scale 'linear' (default). 'log' if input x is log10 values of x, to compute on small x values. + % + % Outputs: + % f computed pdf. + % (if ray method with vpa) can be symbolic for small values + % (if ellipse method with 'log' x_scale) log10 of f + % f_err (if ray method with Monte-Carlo integration) standard error of the output f + % (if Ruben's method) upper error bound of the output f + % (if Imhof's method) logical array, true for outputs too close to 0 or 1 to compute exactly with + % default settings. Try stricter tolerances. + % (if ellipse method) relative error of f + % (if ellipse method with 'log' x_scale) log10 of relative error of f + % x_grid if input x is 'full', this returns the x points at which f is returned % % See also: - % Interactive demos - % gx2cdf - % gx2cdf_davies - % gx2cdf_imhof - % gx2cdf_ruben + % Getting Started guide parser = inputParser; + parser.KeepUnmatched=true; addRequired(parser,'x',@(x) isreal(x) || strcmpi(x,'full')); addRequired(parser,'w',@(x) isreal(x) && isrow(x)); addRequired(parser,'k',@(x) isreal(x) && isrow(x)); addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); - addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); addRequired(parser,'s',@(x) isreal(x) && isscalar(x)); + addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); addParameter(parser,'method','auto'); 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,'xrange',[],@(x) isscalar(x) && (x>0)); - [~,v]=gx2stat(w,k,lambda,m,s); + [~,v]=gx2stat(w,k,lambda,s,m); addParameter(parser,'n_grid',1e4+1,@(x) isscalar(x) && (x>0)); + addParameter(parser,'diff',false,@islogical); addParameter(parser,'dx',sqrt(v)/1e4,@(x) isreal(x) && isscalar(x) && (x>=0)); % default derivative step-size is sd/100. - parse(parser,x,w,k,lambda,m,s,varargin{:}); + parse(parser,x,w,k,lambda,s,m,varargin{:}); method=parser.Results.method; + diff_flag=parser.Results.diff; if strcmpi(x,'full') - method='conv'; + method='ifft'; end - if ~s && length(unique(w))==1 && ~strcmpi(x,'full') - f=ncx2pdf((x-m)/unique(w),sum(k),sum(lambda))/abs(unique(w)); - elseif strcmpi(method,'ifft') - [f,xgrid]=gx2pdf_ifft(x,w,k,lambda,m,s,varargin{:}); - elseif strcmpi(method,'conv') - [f,xgrid]=gx2pdf_conv(x,w,k,lambda,m,s,varargin{:}); - elseif strcmpi(method,'diff') - p_left=gx2cdf(x-dx,w,k,lambda,m,s,varargin{:}); - p_right=gx2cdf(x+dx,w,k,lambda,m,s,varargin{:}); + if ~diff_flag + if strcmpi(method,'auto') + if ~s && length(unique(w))==1 && ~strcmpi(x,'full') + f=ncx2pdf((x-m)/unique(w),sum(k),sum(lambda))/abs(unique(w)); + elseif sum(w)==0 && s % only normal term + f=normpdf(x,m,s); + else + f=gx2_imhof(x,w,k,lambda,s,m,varargin{:},'output','pdf'); + end + elseif strcmpi(method,'imhof') + f=gx2_imhof(x,w,k,lambda,s,m,varargin{:},'output','pdf'); + elseif strcmpi(method,'ruben') + if s || ~(all(w>0)||all(w<0)) + error("Ruben's method can only be used when all w are the same sign and s=0.") + else + f=gx2_ruben(x,w,k,lambda,m,varargin{:},'output','pdf'); + end + elseif strcmpi(method,'ellipse') + if s || ~(all(w>0)||all(w<0)) + error("The ellipse approximation can only be used when all w are the same sign and s=0.") + else + [f,f_err]=gx2_ellipse(x,w,k,lambda,m,varargin{:},'output','pdf'); + end + elseif strcmpi(method,'ray') + [f,f_err]=gx2pdf_ray(x,w,k,lambda,s,m,varargin{:}); + elseif strcmpi(method,'ifft') + [f,xgrid]=gx2_ifft(x,w,k,lambda,s,m,varargin{:},'output','pdf'); + f_err=[]; + end + else % if numerical differentiation + dx=parser.Results.dx; + p_left=gx2cdf(x-dx,w,k,lambda,s,m,varargin{:}); + p_right=gx2cdf(x+dx,w,k,lambda,s,m,varargin{:}); f=(p_right-p_left)/(2*dx); f=max(f,0); end diff --git a/gx2pdf_conv.m b/gx2pdf_conv.m deleted file mode 100644 index d2ebc4d..0000000 --- a/gx2pdf_conv.m +++ /dev/null @@ -1,124 +0,0 @@ -function [f,xgrid]=gx2pdf_conv(x,w,k,lambda,m,s,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,w,k,lambda,m,s) - % f=gx2pdf(x,w,k,lambda,m,s,'dx',1e-3) - % [f,xfull]=gx2pdf('full',w,k,lambda,m,s) - % f=gx2pdf(x,w,k,lambda,m,s,'method','conv','AbsTol',0,'RelTol',1e-7) - % - % Example: - % f=gx2pdf(25,[1 -5 2],[1 2 3],[2 3 7],5,0) - % f=gx2pdf([17 25],[1 -5 2],[1 2 3],[2 3 7],5,0,'method','conv') - % [f,xfull]=gx2pdf('full',[1 -5 2],[1 2 3],[2 3 7],5,0) - % - % Required inputs: - % x points at which to evaluate the pdf, or 'full', to return - % points xfull and f over a full range of the pdf (this uses - % the 'conv' method) - % w row vector of weights of the non-central chi-squares - % k row vector of degrees of freedom of the non-central chi-squares - % lambda row vector of non-centrality paramaters (sum of squares of - % means) of the non-central chi-squares - % m mean of normal term - % s sd of normal term - % - % Optional name-value inputs: - % method 'diff' (default) for differentiating the generalized - % chi-square cdf, or 'conv' for convolving noncentral - % chi-square pdf's - % dx step-size of fineness (for convolving or differentiating) - % AbsTol absolute error tolerance for the output when using 'diff' - % RelTol relative error tolerance for the output when using 'diff' - % The absolute OR the relative tolerance is satisfied. - % - % Output: - % f computed pdf at points, or over full range - % xfull x-values over the full pdf range, returned when x is 'full' - % - % See also: - % Interactive demos - % gx2cdf - % gx2cdf_davies - % gx2cdf_imhof - % gx2cdf_ruben - - parser = inputParser; - addRequired(parser,'x',@(x) isreal(x) || strcmpi(x,'full')); - addRequired(parser,'w',@(x) isreal(x) && isrow(x)); - addRequired(parser,'k',@(x) isreal(x) && isrow(x)); - addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); - addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); - addRequired(parser,'s',@(x) isreal(x) && isscalar(x)); - % range of x over which to compute ifft. - % default span is upto 4 sd from mean - [mu,v]=gx2stat(w,k,lambda,m,s); - addParameter(parser,'xrange',max(abs(mu+[-1 1]*4*sqrt(v))),@(x) isscalar(x) && (x>0)); - % number of grid points for ifft - addParameter(parser,'n_grid',1e4+1,@(x) isscalar(x) && (x>0)); - parse(parser,x,w,k,lambda,m,s,varargin{:}); - - xrange=parser.Results.xrange; - - w_idx=find(w); % non-zero w indices - - if isempty(xrange) % no xrange provided - - % find the xrange over which to convolve - % xrange is a length that covers the - % widest of the constituent pdfs - xrange=nan(1,nnz(w)); - for i=1:nnz(w) - xrange(i)=gx2inv(1-eps,abs(w(w_idx(i))),k(w_idx(i)),lambda(w_idx(i)),0,0); - end - if s - xrange=[xrange,2*norminv(1-eps,0,s)]; - end - xrange=max(2*xrange); - xrange=xrange-mod(xrange,dx); % to center around 0 - end - - dx=parser.Results.dx; - n_grid=round(parser.Results.n_grid); - % make n_grid odd - if mod(n_grid,2) - n_grid=n_grid+1; - end - xgrid=-xrange:dx:xrange; - xgrid=linspace(-xrange,xrange,n_grid); - dx=2*xrange/n_grid; - - % compute non-central and normal pdf's - % over this span - ncpdfs=nan(nnz(w),length(xgrid)); - for i=1:nnz(w) - pdf=gx2pdf(xgrid,w(w_idx(i)),k(w_idx(i)),lambda(w_idx(i)),0,0); - pdf(isinf(pdf))=max(pdf(~isinf(pdf))); - ncpdfs(i,:)=pdf; - end - if s - ncpdfs=[ncpdfs;normpdf(xgrid,0,abs(s))]; - end - f=ifft(prod(fft(ncpdfs,[],2),1)); - if any(isnan(f))||any(isinf(f)) - error('Convolution method failed. Try differentiation method.') - end - if ~mod(size(ncpdfs,1),2) - f=ifftshift(f); - end - f=f/(sum(f)*dx); - xgrid=xgrid+m; - - if ~strcmpi(x,'full') - F=griddedInterpolant(xgrid,f); - f=F(x); - end -end \ No newline at end of file diff --git a/gx2pdf_ifft.m b/gx2pdf_ifft.m deleted file mode 100644 index ef9bc36..0000000 --- a/gx2pdf_ifft.m +++ /dev/null @@ -1,103 +0,0 @@ -function [f,xgrid]=gx2pdf_ifft(x,w,k,lambda,m,s,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,w,k,lambda,m,s) - % f=gx2pdf(x,w,k,lambda,m,s,'dx',1e-3) - % [f,xfull]=gx2pdf('full',w,k,lambda,m,s) - % f=gx2pdf(x,w,k,lambda,m,s,'method','conv','AbsTol',0,'RelTol',1e-7) - % - % Example: - % f=gx2pdf(25,[1 -5 2],[1 2 3],[2 3 7],5,0) - % f=gx2pdf([17 25],[1 -5 2],[1 2 3],[2 3 7],5,0,'method','conv') - % [f,xfull]=gx2pdf('full',[1 -5 2],[1 2 3],[2 3 7],5,0) - % - % Required inputs: - % x points at which to evaluate the pdf, or 'full', to return - % points xfull and f over a full range of the pdf (this uses - % the 'conv' method) - % w row vector of weights of the non-central chi-squares - % k row vector of degrees of freedom of the non-central chi-squares - % lambda row vector of non-centrality paramaters (sum of squares of - % means) of the non-central chi-squares - % m mean of normal term - % s sd of normal term - % - % Optional name-value inputs: - % method 'diff' (default) for differentiating the generalized - % chi-square cdf, or 'conv' for convolving noncentral - % chi-square pdf's - % dx step-size of fineness (for convolving or differentiating) - % AbsTol absolute error tolerance for the output when using 'diff' - % RelTol relative error tolerance for the output when using 'diff' - % The absolute OR the relative tolerance is satisfied. - % - % Output: - % f computed pdf at points, or over full range - % xfull x-values over the full pdf range, returned when x is 'full' - % - % See also: - % Interactive demos - % gx2cdf - % gx2cdf_davies - % gx2cdf_imhof - % gx2cdf_ruben - - parser = inputParser; - addRequired(parser,'x',@(x) isreal(x) || strcmpi(x,'full')); - addRequired(parser,'w',@(x) isreal(x) && isrow(x)); - addRequired(parser,'k',@(x) isreal(x) && isrow(x)); - addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); - addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); - addRequired(parser,'s',@(x) isreal(x) && isscalar(x)); - % range of x over which to compute ifft. - % default span is upto 4 sd from mean - [mu,v]=gx2stat(w,k,lambda,m,s); - addParameter(parser,'xrange',max(abs(mu+[-1 1]*4*sqrt(v))),@(x) isscalar(x) && (x>0)); - % number of grid points for ifft - addParameter(parser,'n_grid',1e4+1,@(x) isscalar(x) && (x>0)); - parse(parser,x,w,k,lambda,m,s,varargin{:}); - - n_grid=round(parser.Results.n_grid); - % make n_grid odd - if mod(n_grid,2) - n_grid=n_grid+1; - end - - n=(n_grid-1)/2; - idx=-n:n; - - xrange=parser.Results.xrange; - dx=xrange/n; - - dt=2*pi/(n_grid*dx); - -% dt=.01; - t=idx*dt; % grid of points to evaluate characteristic function - - phi=gx2char(-t,w,k,lambda,m,s); - - f=fftshift(ifft(ifftshift(phi)))*n_grid*dt/(2*pi); - - % grid of x over which this pdf is computed -% dx=2*pi/(n_grid*dt); - xgrid=idx*dx; - - % normalize - f=f/(sum(f)*dx); - - if ~strcmpi(x,'full') - F=griddedInterpolant(xgrid,f); - f=F(x); - end - - f=max(f,0); -end \ No newline at end of file diff --git a/gx2pdf_imhof.m b/gx2pdf_imhof.m deleted file mode 100644 index f0ebff5..0000000 --- a/gx2pdf_imhof.m +++ /dev/null @@ -1,85 +0,0 @@ -function [f,flag]=gx2pdf_imhof(x,w,k,lambda,m,s,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, which is an extension of Imhof's [1961] method to include the - % s term, so Imhof's method is not separately available in the toolbox - % any more. - % - % 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,w,k,lambda,m,s) - % p=gx2cdf_davies(x,w,k,lambda,m,s,'upper') - % p=gx2cdf_davies(x,w,k,lambda,m,s,'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 - % w row vector of weights of the non-central chi-squares - % k row vector of degrees of freedom of the non-central chi-squares - % lambda row vector of non-centrality paramaters (sum of squares of - % means) of the non-central chi-squares - % m mean of normal term - % s sd of normal 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 for output(s) 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,'w',@(x) isreal(x) && isrow(x)); - addRequired(parser,'k',@(x) isreal(x) && isrow(x)); - addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); - addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); - addRequired(parser,'s',@(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,w,k,lambda,m,s,varargin{:}); - AbsTol=parser.Results.AbsTol; - RelTol=parser.Results.RelTol; - - % compute the integral - if any(strcmpi(parser.UsingDefaults,'AbsTol')) && any(strcmpi(parser.UsingDefaults,'RelTol')) - imhof_integral=arrayfun(@(x) integral(@(u) gx2pdf_imhof_integrand(u,x,w',k',lambda',m,s),0,inf),x); - f=imhof_integral/(2*pi); - else - syms u - imhof_integral=arrayfun(@(x) vpaintegral(@(u) gx2pdf_imhof_integrand(u,x,w',k',lambda',m,s),... - u,0,inf,'AbsTol',AbsTol,'RelTol',RelTol,'MaxFunctionCalls',inf),x); - f=double(imhof_integral/(2*pi)); - end - - flag = f<0; - if any(flag) - warning('gx2pdf output(s) too close to 0 to compute exactly, so clipping to 0. Check the flag output, and try stricter tolerances.') - f=max(f,0); - end - - -end \ No newline at end of file diff --git a/gx2pdf_imhof_integrand.m b/gx2pdf_imhof_integrand.m deleted file mode 100644 index 1aa288f..0000000 --- a/gx2pdf_imhof_integrand.m +++ /dev/null @@ -1,6 +0,0 @@ -function f=gx2pdf_imhof_integrand(u,x,w,k,lambda,m,s) - % define the Imhof integrand for the pdf (w, k, lambda must be column vectors here) - theta=sum(k.*atan(w*u)+(lambda.*(w*u))./(1+w.^2*u.^2),1)/2+u*(m-x)/2; - rho=prod(((1+w.^2*u.^2).^(k/4)).*exp(((w.^2*u.^2).*lambda)./(2*(1+w.^2*u.^2))),1) .* exp(u.^2*s^2/8); - f=cos(theta)./rho; -end \ No newline at end of file diff --git a/gx2pdf_ray.m b/gx2pdf_ray.m new file mode 100644 index 0000000..12fed4c --- /dev/null +++ b/gx2pdf_ray.m @@ -0,0 +1,32 @@ +function [f,f_err]=gx2pdf_ray(x,w,k,lambda,s,m,varargin) + +parser=inputParser; +parser.KeepUnmatched=true; +addRequired(parser,'x',@(x) isreal(x)); +addRequired(parser,'w',@(x) isreal(x) && isrow(x)); +addRequired(parser,'k',@(x) isreal(x) && isrow(x)); +addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); +addRequired(parser,'s',@(x) isreal(x) && isscalar(x)); +addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); +addParameter(parser,'n_rays',1e3); + +parse(parser,x,w,k,lambda,s,m,varargin{:}); + +y=x(:); % flatten input array x into a column vector of levels y + +% find standard quadratic form corresponding to the gx2: +quad=gx2_to_norm_quad_params(w,k,lambda,s,m); +dim=numel(quad.q1); +mu=zeros(dim,1); +v=eye(dim); + +% call int_norm_ray to integrate +[f,f_err]=int_norm_ray(mu,v,quad,'fun_level',y,'output','prob_dens',varargin{:}); + +% reshape flattened array to shape of input x +f=reshape(f,size(x)); +if ~isempty(f_err) + f_err=reshape(f_err,size(x)); +end + +end \ No newline at end of file diff --git a/gx2pdf_ray_fast.m b/gx2pdf_ray_fast.m deleted file mode 100644 index d18c824..0000000 --- a/gx2pdf_ray_fast.m +++ /dev/null @@ -1,98 +0,0 @@ -function f=gx2pdf_ray_fast(x,w,k,lambda,m,s,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, which is an extension of Imhof's [1961] method to include the - % s term, so Imhof's method is not separately available in the toolbox - % any more. - % - % 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,w,k,lambda,m,s) - % p=gx2cdf_davies(x,w,k,lambda,m,s,'upper') - % p=gx2cdf_davies(x,w,k,lambda,m,s,'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 - % w row vector of weights of the non-central chi-squares - % k row vector of degrees of freedom of the non-central chi-squares - % lambda row vector of non-centrality paramaters (sum of squares of - % means) of the non-central chi-squares - % m mean of normal term - % s sd of normal 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 for output(s) 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,'w',@(x) isreal(x) && isrow(x)); - addRequired(parser,'k',@(x) isreal(x) && isrow(x)); - addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); - addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); - addRequired(parser,'s',@(x) isreal(x) && isscalar(x)); - addParameter(parser,'mc_samples',1e4); - - parse(parser,x,w,k,lambda,m,s,varargin{:}); - mc_samples=parser.Results.mc_samples; - - y=x(:); % flatten input array x into a column vector of levels y - - % find standard quadratic form corresponding to the gx2: - quad=gx2_to_norm_quad_params(w,k,lambda,m,s); - dim=numel(quad.q1); - - % uniform random rays (points on n-sphere) - n_z=mvnrnd(zeros(dim,1),eye(dim),mc_samples)'; - n_z=n_z./vecnorm(n_z,2); - - % find the quadratic coefficients across all rays - q2=dot(n_z,quad.q2*n_z); - q1=quad.q1'*n_z; - q0=quad.q0; - - % valid levels y for which there are roots - valid_idx=(sign(q2).*y > sign(q2).*(q0-q1.^2./(4*q2))); - - % discriminant of the quadratic across rays and levels - delta2=q1.^2-4*q2.*(q0-y); % delta^2 - delta=nan(size(delta2)); - delta(valid_idx)=sqrt(delta2(valid_idx)); - - % roots across rays - z=(-q1+cat(3,-1,1).*delta)./(2*q2); - - % 1. phi_ray at all roots across all rays - % 2. sum across two roots on each ray - % 3. divide by quadratic slope - % 4. add across rays, gives us pdf at each level - f=sum(sum(phi_ray(z,dim),3)./delta,2,'omitnan')/mc_samples; - - % reshape flattened array to shape of input x - f=reshape(f,size(x)); -end \ No newline at end of file diff --git a/gx2rnd.m b/gx2rnd.m index d0d0038..5c26ce2 100644 --- a/gx2rnd.m +++ b/gx2rnd.m @@ -1,18 +1,21 @@ -function r=gx2rnd(w,k,lambda,m,s,varargin) +function r=gx2rnd(w,k,lambda,s,m,varargin) - % GX2RND Generates generalized chi-squared random numbers. + % GX2RND Generates generalized chi-squared random numbers. % - % Abhranil Das + % Abhranil Das % Center for Perceptual Systems, University of Texas at Austin + % Comments, questions, bugs to abhranil.das@utexas.edu % If you use this code, please cite: - % A method to integrate and classify normal distributions. + % 1. A method to integrate and classify normal distributions + % 2. New methods for computing the generalized chi-square distribution % % Usage: - % r=gx2rnd(w,k,lambda,m,s) - % r=gx2rnd(w,k,lambda,m,s,sz) - % r=gx2rnd(w,k,lambda,m,s,[sz1,sz2,...]) - % r=gx2rnd(w,k,lambda,m,s,sz,'method',method) + % r=gx2rnd(w,k,lambda,s,m) + % r=gx2rnd(w,k,lambda,s,m,sz) + % r=gx2rnd(w,k,lambda,s,m,[sz1,sz2,...]) + % r=gx2rnd(w,k,lambda,s,m,sz,'method',method) % % Example: % r=gx2rnd([1 -5 2],[1 2 3],[2 3 7],5,1,5) @@ -22,48 +25,45 @@ % k row vector of degrees of freedom of the non-central chi-squares % lambda row vector of non-centrality paramaters (sum of squares of % means) of the non-central chi-squares - % m mean of normal term - % s sd of normal term + % s scale of normal term + % m offset % % Optional inputs: % sz size(s) of the requested array - % method 'sum' (default) to generate non-central chi-square and - % normal numbers and add them. 'norm_quad' to generate - % standard normal vectors and compute their quadratic form. + % method 'sum' (default) to generate non-central chi-square and normal numbers and add them. + % 'norm_quad' to generate standard normal vectors and compute their quadratic form. % % Output: % r random number(s) % % See also: - % Interactive demos - % gx2stat - % gx2_params_norm_quad + % Getting Started guide parser=inputParser; parser.KeepUnmatched=true; addRequired(parser,'w',@(x) isreal(x) && isrow(x)); addRequired(parser,'k',@(x) isreal(x) && isrow(x)); addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); - addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); addRequired(parser,'s',@(x) isreal(x) && isscalar(x)); + addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); addOptional(parser,'sz',@(x) isreal(x)); addParameter(parser,'method','sum',@(s) strcmpi(s,'sum') || strcmpi(s,'norm_quad')); - parse(parser,w,k,lambda,m,s,varargin{:}); + parse(parser,w,k,lambda,s,m,varargin{:}); sz=parser.Results.sz; if isscalar(sz), sz=[sz sz]; end method=parser.Results.method; if strcmpi(method,'sum') - ncxs=arrayfun(@(w,k,lambda) w*ncx2rnd(k,lambda,sz),w,k,lambda,'un',0); + ncxs=arrayfun(@(w,k,lambda) w*ncx2rnd(k,lambda,sz),w,k,lambda,'un',0); - r=zeros(size(ncxs{1})); - for i=1:length(ncxs) - r=r+ncxs{i}; - end - r=r+normrnd(m,s,sz); + r=zeros(size(ncxs{1})); + for i=1:length(ncxs) + r=r+ncxs{i}; + end + r=r+normrnd(m,s,sz); elseif strcmpi(method,'norm_quad') % find the quadratic form of the standard multinormal - quad=gx2_to_norm_quad_params(w,k,lambda,m,s); + quad=gx2_to_norm_quad_params(w,k,lambda,s,m); dim=numel(quad.q1); % generate standard normal vectors diff --git a/gx2stat.m b/gx2stat.m index 1dab5a3..7383b7e 100644 --- a/gx2stat.m +++ b/gx2stat.m @@ -1,12 +1,44 @@ -function [mu,v]=gx2stat(w,k,lambda,m,s) - - parser = inputParser; - addRequired(parser,'w',@(x) isreal(x) && isrow(x)); - addRequired(parser,'k',@(x) isreal(x) && isrow(x)); - addRequired(parser,'lambda',@(x) isreal(x) && isrow(x)); - addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); - addRequired(parser,'s',@(x) isreal(x) && isscalar(x)); - parse(parser,w,k,lambda,m,s); - - mu=dot(w,k+lambda)+m; - v=2*dot(w.^2,k+2*lambda)+s^2; \ No newline at end of file +function [mu,v]=gx2stat(w,k,lambda,s,m) + + % GX2STAT Returns the mean and variance of a generalized chi-squared distribution. + % + % Abhranil Das + % Center for Perceptual Systems, University of Texas at Austin + % Comments, questions, bugs to abhranil.das@utexas.edu + % If you use this code, please cite: + % 1. A method to integrate and classify normal distributions + % 2. New methods for computing the generalized chi-square distribution + % + % Usage: + % [mu,v]=gx2stat(w,k,lambda,s,m) + % + % Example: + % [mu,v]=gx2stat([1 -5 2],[1 2 3],[2 3 7],0,5) + % + % Required inputs: + % w row vector of weights of the non-central chi-squares + % k row vector of degrees of freedom of the non-central chi-squares + % lambda row vector of non-centrality paramaters (sum of squares of + % means) of the non-central chi-squares + % s scale of normal term + % m offset + % + % Outputs: + % mu mean + % v variance + % + % See also: + % Getting Started guide + + parser = inputParser; + addRequired(parser,'w',@(x) isreal(x)); + addRequired(parser,'k',@(x) isreal(x)); + addRequired(parser,'lambda',@(x) isreal(x)); + addRequired(parser,'s',@(x) isreal(x) && isscalar(x)); + addRequired(parser,'m',@(x) isreal(x) && isscalar(x)); + parse(parser,w,k,lambda,s,m); + + mu=dot(w,k+lambda)+m; + v=2*dot(w.^2,k+2*lambda)+s^2; \ No newline at end of file diff --git a/int_norm_ray.m b/int_norm_ray.m new file mode 100644 index 0000000..1c9e1cb --- /dev/null +++ b/int_norm_ray.m @@ -0,0 +1,152 @@ +function [p,p_err,bd_pts]=int_norm_ray(mu,v,dom,varargin) +% Integrate a multinormal distribution over a specified domain, using +% the ray method. + +parser = inputParser; +parser.KeepUnmatched=true; +addRequired(parser,'mu',@isnumeric); +addRequired(parser,'v',@isnumeric); +addRequired(parser,'dom'); +addOptional(parser,'side','lower',@(x) strcmpi(x,'lower') || strcmpi(x,'upper') ); +addParameter(parser,'dom_type','quad'); +addParameter(parser,'output','prob',@(x) strcmpi(x,'prob') || strcmpi(x,'prob_dens') ); +addParameter(parser,'force_mc',false,@islogical); +addParameter(parser,'fun_level',0,@isnumeric); +addParameter(parser,'fun_span',5); +addParameter(parser,'fun_resol',100); +addParameter(parser,'AbsTol',1e-10); +addParameter(parser,'RelTol',1e-2); +addParameter(parser,'vpa',false,@islogical); +addParameter(parser,'n_rays',500); +addParameter(parser,'add_bd_pts',false); +addParameter(parser,'gpu_batch',4e7); + +parse(parser,mu,v,dom,varargin{:}); +output=parser.Results.output; +fun_level=parser.Results.fun_level; +force_mc=parser.Results.force_mc; +AbsTol=parser.Results.AbsTol; +RelTol=parser.Results.RelTol; +n_rays=parser.Results.n_rays; +vpaflag=parser.Results.vpa; +dom_type=parser.Results.dom_type; +gpu_batch=parser.Results.gpu_batch; +isgpu=canUseGPU(); +dim=length(mu); + +if parser.Results.add_bd_pts + global bd_pts +end +bd_pts=[]; + +if force_mc||dim>4 + % Monte-Carlo integration + + if strcmpi(dom_type,'quad')&&~vpaflag&&isgpu&&gpu_batch % batch-process on GPU + disp('Using GPU. If this is slower, set gpu_batch to 0.'); + gpudev=gpuDevice; + reset(gpudev) + + gpu_batch_adj=round(gpu_batch/numel(fun_level)); % adjusted GPU batch size given no. of simultaneous evaluation points + n_batches=ceil(n_rays/gpu_batch_adj); % no. of gpu batches + n_rays_batches=gpu_batch_adj*ones(1,n_batches); % no. of rays in each batch + remainder=rem(n_rays,gpu_batch_adj); + if remainder, n_rays_batches(end)=remainder; end + p_batches=nan(numel(fun_level),n_batches); + p2_batches=nan(numel(fun_level),n_batches); % p squared (for SEM calculation) + for i=1:n_batches + n_z=mvnrnd(zeros(dim,1),eye(dim),n_rays_batches(i))'; + try + n_z=gpuArray(n_z); + p_rays=norm_prob_across_rays(mu,v,dom,n_z,varargin{:}); + p_batches(:,i)=gather(sum(p_rays,2)); + p2_batches(:,i)=gather(sum(p_rays.^2,2)); + catch errmsg + warning('GPU error. If GPU is out of memory, try setting gpu_batch to a smaller value, or to 0 to not use GPU. See error below.') + rethrow(errmsg) + end + end + % mean across all batches + p=sum(p_batches,2)/n_rays; + % mean square across all batches + p2=sum(p2_batches,2)/n_rays; + % SEM of p + p_err=sqrt((p2-p.^2)/n_rays); + + else % process all on CPU + % uniform random rays (points on n-sphere) + n_z=mvnrnd(zeros(dim,1),eye(dim),n_rays)'; + if ~vpaflag + p_rays=norm_prob_across_rays(mu,v,dom,n_z,varargin{:}); + p=mean(p_rays,2); + p_err=std(p_rays,0,2)/sqrt(n_rays); + else % use VPA + [p_rays,~,p_sym_sum]=norm_prob_across_rays(mu,v,dom,n_z,varargin{:}); + p_sum=sum(p_rays,2); + + % probabilities too small even with sym: + sym_idx=cellfun(@(x) isa(x,'sym'),p_sym_sum); + symfail_idx=logical(vpa(p_sym_sum(sym_idx))==0); + if any(symfail_idx) + warning('Probability on some rays too small even for variable precision. Returning 0.') + sym_indices=find(sym_idx); + p_sym_sum(sym_indices(symfail_idx))={0}; + end + + % include the symbolic only if it's + % bigger than RelTol of the numeric + log_p_sym_sum=cellfun(@(x) double(log10(x)),p_sym_sum); + sym_shortlist_idx= sym_idx & (log_p_sym_sum > log10(p_sum) + log10(RelTol)); + if any(sym_shortlist_idx) + warning('Probability too small for double precision. Returning as symbol, use vpa to evaluate.') + % create cell array, combining numeric and symbolic + % probabilitites + p=num2cell(p_sum); + % p(sym_shortlist_idx)=num2cell(p_sum(sym_shortlist_idx)+cell2sym(p_sym_sum(sym_shortlist_idx))); + p(sym_shortlist_idx)=num2cell(cell2sym(p_sym_sum(sym_shortlist_idx))); + p=cellfun(@(x) x/n_rays,p,'un',0); + % p=(p_sum+p_sym_sum)/numel(p_rays); + p_err=[]; + else + p=mean(p_rays,2); + p_err=std(p_rays,0,2)/sqrt(n_rays); + end + end + end + + if strcmpi(output,'prob') + if isnumeric(p) + p=p/2; + elseif iscell(p) + p=cellfun(@(x) x/2,p,'un',0); + end + end + +else + + % grid integration + fun_level=parser.Results.fun_level; + + if dim==1 + p=norm_prob_across_angles(mu,v,dom,varargin{:}); + elseif dim==2 + if numel(fun_level)==1 % if integral needed only at one function level, + % integrate with 'arrayvalued', false to evaluate + % simultaneously across grid of theta and speed up + p=integral(@(theta) norm_prob_across_angles(mu,v,dom,varargin{:},'theta',theta),0,pi,'AbsTol',AbsTol,'RelTol',RelTol); + else + % integrate with 'arrayvalued', true to simultaneously + % integrate at multiple function levels + p=integral(@(theta) norm_prob_across_angles(mu,v,dom,varargin{:},'theta',theta),0,pi,'arrayvalued',true,'AbsTol',AbsTol,'RelTol',RelTol); + end + elseif dim==3 + % use arrayfun to integrate at each fun_level separately, + % because integral2 cannot integrate vector-valued function. + p=arrayfun(@(f) integral2(@(theta,phi) norm_prob_across_angles(mu,v,dom,varargin{:},'theta',theta,'phi',phi,'fun_level',f),0,pi/2,0,2*pi,'AbsTol',AbsTol,'RelTol',RelTol),fun_level); + elseif dim==4 + % use arrayfun to integrate at each fun_level separately, + % because integral3 cannot integrate vector-valued function. + p=arrayfun(@(f) integral3(@(theta,phi,psi) norm_prob_across_angles(mu,v,dom,varargin{:},'theta',theta,'phi',phi,'psi',psi,'fun_level',f),0,pi/2,0,2*pi,0,pi,'AbsTol',AbsTol,'RelTol',RelTol),fun_level); + end + p_err=[]; +end \ No newline at end of file diff --git a/log_gx2cdf.m b/log_gx2cdf.m new file mode 100644 index 0000000..e9bb22c --- /dev/null +++ b/log_gx2cdf.m @@ -0,0 +1,13 @@ +function l=log_gx2cdf(x,w,k,lambda,s,m,varargin) +p=gx2cdf(x,w,k,lambda,s,m,varargin{:}); +if isnumeric(p) + if p<=0 + l=p; + elseif p>0 + l=log10(p); + end +elseif iscell(p) + l=double(log10(cell2sym(p))); +end +% end +end \ No newline at end of file diff --git a/norm_prob_across_angles.m b/norm_prob_across_angles.m new file mode 100644 index 0000000..fbd9179 --- /dev/null +++ b/norm_prob_across_angles.m @@ -0,0 +1,75 @@ +function [p_angles,bd_pts_angles]=norm_prob_across_angles(mu,v,dom,varargin) + + % parse inputs + parser=inputParser; + parser.KeepUnmatched=true; + addRequired(parser,'mu',@isnumeric); + addRequired(parser,'v',@isnumeric); + addRequired(parser,'dom',@(x) isstruct(x) || isa(x,'function_handle') || ismatrix(x)); + addOptional(parser,'side','lower',@(x) strcmpi(x,'lower') || strcmpi(x,'upper') ); + addParameter(parser,'dom_type','quad'); + addParameter(parser,'output','prob'); % probability or probability density + addParameter(parser,'fun_span',5); + addParameter(parser,'fun_resol',100); + addParameter(parser,'n_bd_pts',500); + addParameter(parser,'add_bd_pts',false); + addParameter(parser,'theta',nan,@isnumeric); + addParameter(parser,'phi',nan,@isnumeric); + addParameter(parser,'psi',nan,@isnumeric); + + parse(parser,mu,v,dom,varargin{:}); + + n_bd_pts=parser.Results.n_bd_pts; + theta=parser.Results.theta; + phi=parser.Results.phi; + psi=parser.Results.psi; + output=parser.Results.output; + + dim=length(mu); + if dim>1 && any(strcmpi(parser.UsingDefaults,'theta')) && any(strcmpi(parser.UsingDefaults,'phi')) + points=fibonacci_sphere(n_bd_pts); + [az,el]=cart2sph(points(1,:),points(2,:),points(3,:)); + theta=pi/2-el; + phi=az; + if dim==2 + theta=linspace(0,pi,n_bd_pts); + end + end + + if dim==1 + n_z=1; + elseif dim==2 + [x,y]=pol2cart(theta,1); + n_z=[x;y]; + elseif dim==3 + az=phi; + el=pi/2-theta; + [x,y,z]=sph2cart(az,el,1); + n_z=[x(:)';y(:)';z(:)']; + elseif dim==4 + x=cos(psi); + y=sin(psi).*cos(theta); + z=sin(psi).*sin(theta).*cos(phi); + w=sin(psi).*sin(theta).*sin(phi); + n_z=[x(:)';y(:)';z(:)';w(:)']; + end + + [p_angles,bd_pts_angles]=norm_prob_across_rays(mu,v,dom,n_z,varargin{:}); + + % factors: first multiply by 2/total angle + if dim==2 + p_angles=p_angles/pi; + elseif dim==3 + p_angles=reshape(p_angles,size(theta)); + p_angles=p_angles.*sin(theta)/(2*pi); + elseif dim==4 + p_angles=reshape(p_angles,size(theta)); + p_angles=p_angles.*sin(theta).*sin(psi).^2/(pi^2); + end + + if strcmpi(output,'prob') + % for prob. density, the factor is 2/total angle, + % but for prob, it's 1/total angle. + % so extra factor of 1/2. + p_angles=p_angles/2; + end \ No newline at end of file diff --git a/norm_prob_across_rays.m b/norm_prob_across_rays.m new file mode 100644 index 0000000..94a233b --- /dev/null +++ b/norm_prob_across_rays.m @@ -0,0 +1,93 @@ +function [p_rays,bd_pts_rays,p_sym_sum,sym_idx]=norm_prob_across_rays(mu,v,dom,n_z,varargin) + +% parse inputs +parser=inputParser; +parser.KeepUnmatched=true; +addRequired(parser,'mu',@isnumeric); +addRequired(parser,'v',@isnumeric); +addRequired(parser,'dom',@(x) isstruct(x) || isa(x,'function_handle') || ismatrix(x)); +addRequired(parser,'n_z',@isnumeric); +addOptional(parser,'side','lower',@(x) strcmpi(x,'lower') || strcmpi(x,'upper') ); +addParameter(parser,'output','prob',@(x) strcmpi(x,'prob') || strcmpi(x,'prob_dens') ); +addParameter(parser,'dom_type','quad'); +addParameter(parser,'vpa',false,@islogical); +addParameter(parser,'fun_level',0,@isnumeric); +addParameter(parser,'fun_span',5); +addParameter(parser,'fun_resol',100); +addParameter(parser,'fun_grad',[],@(x) isa(x,'function_handle')); +addParameter(parser,'n_bd_pts',1e4); +addParameter(parser,'add_bd_pts',false); + +parse(parser,mu,v,dom,n_z,varargin{:}); + +dom_type=parser.Results.dom_type; +fun_level=parser.Results.fun_level; +output=parser.Results.output; +vpaflag=parser.Results.vpa; +dim=length(mu); + +n_z=n_z./vecnorm(n_z,2,1); + +% ray-trace if necessary +if ~strcmpi(dom_type,'quad') || parser.Results.add_bd_pts + % ray-trace the standardized domain/function + dom_standard_raytrace=@(n) standard_ray_trace(dom,n,varargin{:},'mu',mu,'v',v); + + % initial signs and boundary distances in standardized space + [init_sign,z]=dom_standard_raytrace(n_z); +end + +if strcmpi(dom_type,'quad') + % standardized boundary coefficients + quad_s=standard_quad(dom,mu,v); + + if ~vpaflag + p_rays=gx2_ray_integrand(fun_level,n_z,quad_s,varargin{:}); + else + [p_rays,p_sym_sum,sym_idx]=gx2_ray_integrand(fun_level,n_z,quad_s,varargin{:}); + end +else + if strcmpi(output,'prob') + % probability on rays + if ~vpaflag + p_rays=cellfun(@(init_sign_ray,z_ray) prob_ray(init_sign_ray,z_ray,dim,varargin{:}),num2cell(init_sign),z); + % if there are roots on rays with 0 prob, + % notify to turn on vpa + n_roots=cellfun(@(z) numel(z)>0,z); + if nnz(n_roots&(~p_rays)) + warning("Some rays contain probabilities smaller than realmin=1e-308, returning 0. Set 'vpa' to true to compute these with variable precision.") + end + else + p_rays=cellfun(@(init_sign_ray,z_ray) prob_ray(init_sign_ray,z_ray,dim,varargin{:}),num2cell(init_sign),z,'un',0); + end + elseif strcmpi(output,'prob_dens') % probability density calculations + + % gradient of standardized function + fun_grad=parser.Results.fun_grad; + standard_gradf=@(z) standard_fun_grad(z,fun_grad,mu,v); + + % probability density on rays + p_rays=cellfun(@(n_ray,z_ray) prob_dens_ray(n_ray,z_ray,standard_gradf), num2cell(n_z,1),z); + end + + % if p_rays is a cell but there are no symbols, convert to numeric array + if iscell(p_rays) + num_idx=cellfun(@isnumeric, p_rays); + if all(num_idx) + p_rays=cell2mat(p_rays); + end + end +end + +if parser.Results.add_bd_pts + % standard boundary points + std_bd_pts_ray=cellfun(@(z_ray,n_ray) z_ray.*n_ray, z,num2cell(n_z,1),'un',0); + + % boundary points + bd_pts_rays=sqrtm(v)*horzcat(std_bd_pts_ray{:})+mu; + + global bd_pts + bd_pts=[bd_pts,bd_pts_rays]; +else + bd_pts_rays=[]; +end diff --git a/norm_quad_to_gx2_params.m b/norm_quad_to_gx2_params.m index 570151a..c225863 100644 --- a/norm_quad_to_gx2_params.m +++ b/norm_quad_to_gx2_params.m @@ -1,15 +1,18 @@ -function [w,k,lambda,m,s]=norm_quad_to_gx2_params(mu,v,quad) +function [w,k,lambda,s,m]=norm_quad_to_gx2_params(mu,v,quad) % NORM_QUAD_TO_GX2_PARAMS A quadratic form of a normal vector is distributed % as a generalized chi-squared. This function takes the multinormal parameters - % and the quadratic coeffs and returns the parameters of the generalized + % and the quadratic coefficients 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. + % Abhranil Das + % Center for Perceptual Systems, University of Texas at Austin + % Comments, questions, bugs to abhranil.das@utexas.edu + % If you use this code, please cite: + % 1. A method to integrate and classify normal distributions + % 2. New methods for computing the generalized chi-square distribution % % Example: % mu=[1;2]; % mean @@ -19,7 +22,7 @@ % quad.q1=[-1;0]; % quad.q0=-1; % - % [w,k,lambda,m,s]=gx2_params_norm_quad(mu,v,quad) + % [w,k,lambda,s,m]=gx2_params_norm_quad(mu,v,quad) % % Required inputs: % mu column vector of normal mean @@ -34,15 +37,11 @@ % k row vector of degrees of freedom of the non-central chi-squares % lambda row vector of non-centrality paramaters (sum of squares of % means) of the non-central chi-squares - % m mean of normal term - % s sd of normal term + % s scale of normal term + % m offset % % See also: - % Interactive demos - % gx2rnd - % gx2stat - % gx2cdf - % gx2pdf + % Getting Started guide % standardize the space q2s=0.5*(quad.q2+quad.q2'); % symmetrize q2 diff --git a/paper 2/1_cdf.pdf b/paper 2/1_cdf.pdf deleted file mode 100644 index 988a923..0000000 Binary files a/paper 2/1_cdf.pdf and /dev/null differ diff --git a/paper 2/1a_gx2cdf.pdf b/paper 2/1a_gx2cdf.pdf deleted file mode 100644 index 31d83c3..0000000 Binary files a/paper 2/1a_gx2cdf.pdf and /dev/null differ diff --git a/paper 2/1a_gx2pdf.pdf b/paper 2/1a_gx2pdf.pdf deleted file mode 100644 index acf24f5..0000000 Binary files a/paper 2/1a_gx2pdf.pdf and /dev/null differ diff --git a/paper 2/1a_ncx2_compare.fig b/paper 2/1a_ncx2_compare.fig deleted file mode 100644 index 0189187..0000000 Binary files a/paper 2/1a_ncx2_compare.fig and /dev/null differ diff --git a/paper 2/1a_norm_quad.pdf b/paper 2/1a_norm_quad.pdf deleted file mode 100644 index 44e4129..0000000 Binary files a/paper 2/1a_norm_quad.pdf and /dev/null differ diff --git a/paper 2/1b_ncx2_compare.pdf b/paper 2/1b_ncx2_compare.pdf deleted file mode 100644 index c58bb7d..0000000 Binary files a/paper 2/1b_ncx2_compare.pdf and /dev/null differ diff --git a/paper 2/1c_gx2_compare.fig b/paper 2/1c_gx2_compare.fig deleted file mode 100644 index cb49669..0000000 Binary files a/paper 2/1c_gx2_compare.fig and /dev/null differ diff --git a/paper 2/1c_gx2_compare.pdf b/paper 2/1c_gx2_compare.pdf deleted file mode 100644 index 3632001..0000000 Binary files a/paper 2/1c_gx2_compare.pdf and /dev/null differ diff --git a/paper 2/1c_imhof_integrand.pdf b/paper 2/1c_imhof_integrand.pdf deleted file mode 100644 index 0cb5814..0000000 Binary files a/paper 2/1c_imhof_integrand.pdf and /dev/null differ diff --git a/paper 2/2_gx2pdf.pdf b/paper 2/2_gx2pdf.pdf deleted file mode 100644 index 01641af..0000000 Binary files a/paper 2/2_gx2pdf.pdf and /dev/null differ diff --git a/paper 2/2a_davies_pdf.fig b/paper 2/2a_davies_pdf.fig deleted file mode 100644 index d8ee6ed..0000000 Binary files a/paper 2/2a_davies_pdf.fig and /dev/null differ diff --git a/paper 2/2a_ray_pdf_cartoon.pdf b/paper 2/2a_ray_pdf_cartoon.pdf deleted file mode 100644 index 90bb703..0000000 Binary files a/paper 2/2a_ray_pdf_cartoon.pdf and /dev/null differ diff --git a/paper 2/2b_polyfun_pdf.pdf b/paper 2/2b_polyfun_pdf.pdf deleted file mode 100644 index 153b762..0000000 Binary files a/paper 2/2b_polyfun_pdf.pdf and /dev/null differ diff --git a/paper 2/2c_gx2pdf_compare.pdf b/paper 2/2c_gx2pdf_compare.pdf deleted file mode 100644 index c2f1950..0000000 Binary files a/paper 2/2c_gx2pdf_compare.pdf and /dev/null differ diff --git a/paper 2/3_dprime_check.fig b/paper 2/3_dprime_check.fig deleted file mode 100644 index f437980..0000000 Binary files a/paper 2/3_dprime_check.fig and /dev/null differ diff --git a/paper 2/3_dprime_check.pdf b/paper 2/3_dprime_check.pdf deleted file mode 100644 index d142140..0000000 Binary files a/paper 2/3_dprime_check.pdf and /dev/null differ diff --git a/paper 2/dprime_check.mat b/paper 2/dprime_check.mat deleted file mode 100644 index a89c7ad..0000000 Binary files a/paper 2/dprime_check.mat and /dev/null differ diff --git a/paper 2/plots.m b/paper 2/plots.m deleted file mode 100644 index d596475..0000000 --- a/paper 2/plots.m +++ /dev/null @@ -1,471 +0,0 @@ -%% define colours -orange=[255 123 49]/255; -sky=[0 179 255]/255; -cobalt=[0 116 255]/255; - -%% quadratic picture corresponding to gx2cdf - -w=[1 -1]; -k=[1 1]; -lambda=[2 4]; -m=0; -s=0; - -x_full=[linspace(-20,-.5,200) linspace(-.5,.5,100) linspace(.5,20,50)]; -x=linspace(-15,15,7); - -colors=spring(length(x)); - -% plot pdf -f_full=gx2pdf(x_full,w,k,lambda,m,s,'method','conv'); -figure; hold on - -int_idx=x_full1e-3; - -f_ray_plot=f_ray; -f_ray_plot(mid_idx)=f_ray_plot(mid_idx)*1e4-13; -f_ray_plot(~mid_idx)=log10(f_ray_plot(~mid_idx)); - -f_conv_plot=f_conv; -f_conv_plot(mid_idx)=f_conv_plot(mid_idx)*1e4-13; -f_conv_plot(~mid_idx)=log10(f_conv_plot(~mid_idx)); - -f_im_plot=f_im; -f_im_plot(mid_idx)=f_im_plot(mid_idx)*1e4-13; -f_im_plot(~mid_idx)=log10(f_im_plot(~mid_idx)); - -figure; hold on -plot(x_plot,f_ray_plot,'-o','markersize',6,'color',sky,'markerfacecolor','w','linewidth',1) -plot(x_plot,f_conv_plot,'ok','marker','.','markersize',5,'linewidth',1) -plot(x_plot,f_im_plot,'or','markersize',3,'color',orange,'linewidth',1) - -yline([-3 log10(realmin)]) -x_trans_plot=x_plot(find(diff(mid_idx))); -xline(x_trans_plot) -x_top_ticks=x_mid_plot(end)+20*(log10([1e3 1e4])-log10(x_top(1))); -xlim([-80 x_top_ticks(end)]) -set(gca,'xtick',[-80 -60 x_trans_plot x_top_ticks]) -set(gca,'XTickLabel',{'-10^4','-10^3','-47','22','10^3','10^4'}) -ylim([-400 .04*1e4-13]) -set(gca,'ytick',[-400 -200 -3 .04*1e4-13]) -set(gca,'YTickLabel',{'10^{-400}','10^{-200}','0.001','0.04'}) - -xlabel("$\tilde{\chi}$",'interpreter','latex') -ylabel 'pdf' -set(gca,'fontsize',13) - -%% extended d' check with vpa - -mu_1=[0;0;0]; - -% both normals have the same covariance, so we can check against -% the true d' (Mahalanobis distance). -v=[1 .5 .7; - .5 2 1 ; - .7 1 3]; - -steps=logspace(0,5,20); -d_true=nan(size(steps)); -d_gx2=nan(size(steps)); -d_ray=nan(size(steps)); - -d_realmin=-2*norminv(realmin); - -for i=1:length(steps) - i - mu_2=steps(i)*[1;1;1]; - - results_gx2=classify_normals([mu_1,v],[mu_2,v],'method','gx2','plotmode',false); - d_true(i)=results_gx2.norm_d_e; - d_gx2(i)=results_gx2.norm_d_b; - - tic - if d_true(i)d_realmin,1); -plot(d_true(first_vpa-1:end),rel_err_ray(first_vpa-1:end),'-o','color',cobalt,'markersize',4,'markerfacecolor',cobalt) -plot(d_true(1:first_vpa-1),rel_err_ray(1:first_vpa-1),'-o','color',sky,'markersize',4,'markerfacecolor',sky) -set(gca,'xscale','log') -set(gca,'yscale','log') -ylabel("rel. error $ | \hat{d}' - d' | / d' $",'Interpreter','latex'); -xline([d_realmin d_symmin]) % largest computable d', corr. to the smallest possible error representable in double-precision -% yline(eps) % machine epsilon for double precision -xlabel("true $d'$",'interpreter','latex') -ylim([1e-18 1e-2]) -set(gca,'ytick',[1e-18 1e-15 1e-10 1e-5]) -set(gca,'yticklabel',{'0','10^{-15}','10^{-10}','10^{-5}'}) -set(gca,'fontsize',13) \ No newline at end of file diff --git a/phi_ray.m b/phi_ray.m new file mode 100644 index 0000000..ef8ef72 --- /dev/null +++ b/phi_ray.m @@ -0,0 +1,6 @@ +function f=phi_ray(z,dim) + % standard multinormal density function along a ray + % f=abs(z).*chi2pdf(z.^2,dim); + z2=z.^2; + f=abs(z).*(z2.^(dim/2-1)).*exp(-z2/2)/(2^(dim/2)*gamma(dim/2)); +end \ No newline at end of file diff --git a/standard_quad.m b/standard_quad.m new file mode 100644 index 0000000..2c3d7a6 --- /dev/null +++ b/standard_quad.m @@ -0,0 +1,5 @@ +function quad_s=standard_quad(quad,mu,v) +% standardize quadratic coefficients +quad_s.q2=sqrtm(v)*quad.q2*sqrtm(v); +quad_s.q1=sqrtm(v)*(2*quad.q2*mu+quad.q1); +quad_s.q0=mu'*quad.q2*mu+quad.q1'*mu+quad.q0; \ No newline at end of file diff --git a/test.m b/test.m deleted file mode 100644 index bb23f49..0000000 --- a/test.m +++ /dev/null @@ -1,135 +0,0 @@ - -%% debug IntClassNorm -mu_1=[4; 5]; v_1=[2 1; 1 1]; -mu_2=mu_1; v_2=3*[2 -1; -1 1]; - -results_gx2=classify_normals([mu_1,v_1],[mu_2,v_2],'method','gx2') - -%% pdf of vector function -mu=[-1; -1]; v=[1 0.5; 0.5 2]; -fun=@(x,y) x.*sin(y) - y.*cos(x); - -fun_grad=@(x) [sin(x(2))+x(2)*sin(x(1)); x(1)*cos(x(2))-cos(x(1))]; - -x=linspace(-20,20,50); -% p=norm_fun_cdf(x,mu,v,fun); -% plot(x,p) -tic -f_diff=norm_fun_pdf(x,mu,v,fun,'pdf_method','diff','dx',0.1,'AbsTol',0,'RelTol',1e-5); -time_diff=toc - -tic -f_ray=norm_fun_pdf(x,mu,v,fun,'pdf_method','ray','fun_grad',fun_grad,'AbsTol',0,'RelTol',1e-5); -time_ray=toc -figure; hold on -plot(x,f_diff) -plot(x,f_ray) - -%% develop gx2pdf_ray - -% ncx2 -w=1; -k=4; -lambda=5; -m=0; -s=0; - -x=linspace(0,50,1e3); -f_ray=gx2pdf_ray(x,w,k,lambda,m,s); - -f_nc=ncx2pdf(x,k,lambda); - -figure; hold on -plot(x,f_ray,'-k') -plot(x,f_nc) - -% gx2 - -w=[1 -1]; -k=[1 1]; -lambda=[2 4]; -m=0; -s=0; - -x=linspace(-50,50,1e2); -tic -f_ray=gx2pdf_ray(x,w,k,lambda,m,s); -t_ray=toc; - -tic -f_im=gx2pdf_imhof(x,w,k,lambda,m,s); -t_im=toc; - -figure; hold on -plot(x,f_ray,'-k') -plot(x,f_im) - - -%% IFFT of characteristic function -% w=[-5 1 2]; -% k=[1 2 1]; -% lambda=[0 0 0]; -% m=0; -% s=0; - -w=1; -k=1; -lambda=0; -m=0; -s=0; - -% w=0; -% k=1; -% lambda=0; -% m=0; -% s=1; - -[f_conv,x_conv]=gx2pdf('full',w,k,lambda,m,s,'method','conv','xrange',150,'n_grid',1e6); -[f_ifft,x_ifft]=gx2pdf_ifft('full',w,k,lambda,m,s,'xrange',150,'n_grid',1e6); - -figure; -hold on -plot(x_ifft,f_ifft) -plot(x_conv,f_conv) -set(gca,'yscale','log') -% xlim([-300 300]) -% ylim([1e-20 1]) - -%% vpa for ray integration up to 3D - -format long -dprime_true=100 - -mu_1=[0;0]; v=eye(2); -mu_2=dprime_true*[1;0]; - -results=classify_normals([mu_1,v],[mu_2,v],'AbsTol',0,'RelTol',1e-2,'vpa',true) - -%% sym roots for SO -theta=-1; % theta = -1/0/1 give 2/1/0 valid roots respectively - -q2=2*theta; -q1=theta+1; -q0=1; - -r=valid_quad_roots(q2,q1,q0) - -syms theta - -q2=2*theta; -q1=theta+1; -q0=1; - -r=valid_quad_roots(q2,q1,q0) - -%% simplify psi(n) - -format long -dprime_true=70 - -mu_1=[0;0;0]; v=eye(3); -mu_2=dprime_true*[1;0;0]; - -results=classify_normals([mu_1,v],[mu_2,v],'AbsTol',0,'RelTol',1e-2,'plotmode',0) -dprime=results.norm_d_b -