Updates

Phibar_ray_split.m

@@ -1,19 +1,20 @@
-function [f_big,f_small]=Phibar_ray_split(z,dim)
+function [Phibar_big,Phibar_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); % 2 when z<=-z_c, 1 when -z_c<z<z_c, 0 when z>=z_c
+    z_c=sqrt(chi2inv(0.5,dim)); % if z^2> z_c^2, we'll use chi2cdf('upper').
+    Phibar_big=(z<=-z_c)+(z<z_c); % 2 when z<=-z_c, 1 when -z_c<z<z_c, 0 when 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)&(z<z_c)))); % if z<=-z_c or 0<z<z_c, invert sign
+    Phibar_small=gammainc(z.^2/2,dim/2); %chi2cdf(z.^2,dim);
+    Phibar_small_upper=gammainc(z.^2/2,dim/2,'upper'); %chi2cdf(z.^2,dim,'upper');
+    Phibar_small(abs(z)>=z_c)=Phibar_small_upper(abs(z)>=z_c);
+
+    Phibar_small=Phibar_small.*(1-2*((z<=-z_c)|((z>0)&(z<z_c)))); % if z<=-z_c or 0<z<z_c, invert sign
 
     % if any z is nan (non-existent root),
     % prob. contribution is 0.
-    f_big(isnan(f_big))=0;
-    f_small(isnan(f_small))=0;
+    Phibar_big(isnan(Phibar_big))=0;
+    Phibar_small(isnan(Phibar_small))=0;
 
 end

Phibar_sym.m

@@ -0,0 +1,11 @@
+function p=Phibar_sym(z,dim)
+% complementary cdf on ray using symbol-friendly igamma
+% first use igamma to construct regularized gamma, using which
+% construct chi distribution cdf, using which construct phibar.
+
+% chi CDF is regularized gamma, built from symbol-friendly igamma
+chiCDF=1-igamma(dim/2,z.^2/2)/gamma(dim/2);
+
+% Phibar from chicdf
+p=(1-sign(z).*chiCDF)/2;
+end

gx2_das.m

@@ -0,0 +1,62 @@
+function [p,errflag]=gx2_das(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,'scale','linear',@(x) strcmpi(x,'linear') || strcmpi(x,'log') );
+
+parse(parser,x,w,k,lambda,s,m,varargin{:});
+side=parser.Results.side;
+
+% first merge into unique w's
+[w,~,ic]=uniquetol(w); % unique non-zero eigenvalues
+k=arrayfun(@(x)sum(k(ic==x)),1:numel(w)); % merged total dof's
+lambda=arrayfun(@(x)sum(lambda(ic==x)),1:numel(w)); % merged total non-centralities
+
+if strcmpi(side,'upper') % upper tail
+    [w_max,max_idx]=max(w.*(w>0));
+else % lower tail
+    [w_max,max_idx]=min(w.*(w<0));
+    % x=abs(x);
+    % sgn=-1;
+end
+sgn=1;
+
+k_max=k(max_idx);
+lambda_max=lambda(max_idx);
+
+w_rest=w([1:max_idx-1, max_idx+1:end]);
+k_rest=k([1:max_idx-1, max_idx+1:end]);
+lambda_rest=lambda([1:max_idx-1, max_idx+1:end]);
+
+a=exp(sgn*m/(2*w_max)+s^2/(8*w_max^2))*...
+    prod(exp((lambda_rest.*w_rest*sgn)./(2*(w_max-sgn*w_rest)))./(1-sgn*w_rest/w_max).^(k_rest/2));
+
+if strcmpi(parser.Results.output,'pdf')
+    if strcmpi(parser.Results.scale,'linear')
+        p=a/abs(w_max)*ncx2pdf(x/w_max,k_max,lambda_max);
+    % elseif strcmpi(parser.Results.scale,'linear')
+        % p=(k_max/2-1)*log10(x)+log10(exp(1))*(sgn*m/(2*w_max)+s^2/(8*w_max^2)-x/(2*w_max))-...
+    % log10(gamma(k_max/2))-(k_max/2)*log10(2*w_max)-sum((k_rest/2).*log10(1-sgn*w_rest/w_max));
+    end
+elseif strcmpi(parser.Results.output,'cdf')
+    marcum=marcumq(sqrt(lambda_max),sqrt(x/w_max),k_max/2);
+    p=a*marcum;
+end
+
+errflag=p<0;
+
+if any(errflag)
+    warning('Tail approximation output(s), i.e. Marcum Q-function values, are too small to correctly compute, so clipping to 0. Check the error flag output.')
+    p=max(p,0);
+end
+
+log10_f=(k_max/2-1)*log10(x)+log10(exp(1))*(sgn*m/(2*w_max)+s^2/(8*w_max^2)-x/(2*w_max))-...
+    log10(gamma(k_max/2))-(k_max/2)*log10(2*w_max)-sum((k_rest/2).*log10(1-sgn*w_rest/w_max));

gx2_ifft.m

@@ -17,7 +17,8 @@
     span=max(abs(mu+[-1 1]*100*sqrt(v)));
 else
     % default span is 100* input span
-    span=100*max(abs(x(:)));
+    % span=100*max(abs(x(:)));
+    span=1e5;
 end
 addParameter(parser,'span',span,@(x) isscalar(x) && (x>0));
 % number of grid points for ifft
@@ -36,7 +37,12 @@
 idx=-n:n;
 dx=span/n;
 
-xgrid=idx*dx; % grid of x over which this cdf/pdf is computed
+if strcmpi(x,'full')
+    x_mid=0;
+else
+    x_mid=(min(x(:))+max(x(:)))/2;
+end
+xgrid=x_mid+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
@@ -64,9 +70,10 @@
     phi=gx2char(-t,w,k,lambda,s,m);
 
     if strcmpi(parser.Results.output,'pdf')
+        phi=phi.*exp(1i*x_mid*dt*idx);
         p=fftshift(ifft(ifftshift(phi)))/dx;
     elseif strcmpi(parser.Results.output,'cdf')
-        phi=phi./(1i*t);
+        phi=phi./(1i*t).*exp(1i*x_mid*dt*idx);
         phi(isinf(phi))=0;
         p=0.5+fftshift(ifft(ifftshift(phi)))/dx;
     end

gx2_ifft_old.m

@@ -0,0 +1,85 @@
+function [p,xgrid]=gx2_ifft_old(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);

gx2_imhof_integrand.m

@@ -1,4 +1,4 @@
-function f=gx2_imhof_integrand(u,x,w,k,lambda,s,m,output)
+function [f,theta,rho]=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);

gx2_pearson.m

@@ -0,0 +1,51 @@
+function p=gx2_pearson(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) && 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') );
+
+parse(parser,x,w,k,lambda,s,m,varargin{:});
+side=parser.Results.side;
+
+% first 3 central moments:
+mu(1)=sum(w.*(k+lambda))+m;
+mu(2)=2*sum(w.^2.*(k+2*lambda))+s^2;
+mu(3)=8*sum(w.^3.*(k+3*lambda));
+
+h=8*mu(2)^3/mu(3)^2;
+
+if mu(3)>0
+    y=(x-mu(1))*sqrt(2*h/mu(2))+h;
+    if strcmpi(parser.Results.output,'cdf')
+        if strcmpi(side,'lower')
+            p=chi2cdf(y,h);
+        elseif strcmpi(side,'upper')
+            p=chi2cdf(y,h,'upper');
+        end
+    elseif strcmpi(parser.Results.output,'pdf')
+        p=sqrt(2*h/mu(2))*chi2pdf(y,h);
+    end
+else
+    mu(1)=-mu(1);
+    mu(3)=-mu(3);
+    x=-x;
+    y=(x-mu(1))*sqrt(2*h/mu(2))+h;
+    if strcmpi(parser.Results.output,'cdf')
+        if strcmpi(side,'lower')
+            p=chi2cdf(y,h,'upper');
+        elseif strcmpi(side,'upper')
+            p=chi2cdf(y,h);
+        end
+    elseif strcmpi(parser.Results.output,'pdf')
+        p=sqrt(2*h/mu(2))*chi2pdf(y,h);
+    end
+end
+
+end