-
Notifications
You must be signed in to change notification settings - Fork 2
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
1 parent
07d3325
commit 3b16ed4
Showing
53 changed files
with
1,276 additions
and
1,533 deletions.
There are no files selected for viewing
Binary file not shown.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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); % 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 | ||
|
|
||
| % if any z is nan (non-existent root), | ||
| % prob. contribution is 0. | ||
| f_big(isnan(f_big))=0; | ||
| f_small(isnan(f_small))=0; | ||
|
|
||
| end |
Binary file not shown.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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)); |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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); |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 |
Oops, something went wrong.