Index: /issm/trunk-jpl/src/m/miscellaneous/normcdf_issm.m =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/normcdf_issm.m (revision 14182) +++ /issm/trunk-jpl/src/m/miscellaneous/normcdf_issm.m (revision 14182) @@ -0,0 +1,9 @@ +% +% wrapper for normcdf to avoid using the matlab statistics toolbox. +% +function [p]=normcdf_issm(x,mu,sigma) + + p=(1.+erf((x-mu)/(sigma*sqrt(2.))))/2.; + +end + Index: /issm/trunk-jpl/src/m/miscellaneous/normfit_issm.m =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/normfit_issm.m (revision 14182) +++ /issm/trunk-jpl/src/m/miscellaneous/normfit_issm.m (revision 14182) @@ -0,0 +1,36 @@ +% +% wrapper for normfit to avoid using the matlab statistics toolbox. +% +function [muhat,sigmahat,muci,sigmaci]=normfit_issm(x,alpha) + +% remove any NaN + + x=x(~isnan(x(:))); + n=length(x); + +% explicitly calculate the moments + + muhat =mean(x); + sigmahat=std(x); + + if (nargout>2) + if ~exist('alpha','var') + alpha=0.05; + end + prob=1.-alpha/2.; + + try + muci(1,1) =muhat-tinv(prob,n-1)*sigmahat/sqrt(n); + muci(2,1) =muhat+tinv(prob,n-1)*sigmahat/sqrt(n); + sigmaci(1,1)=sigmahat*sqrt((n-1)/chi2inv(prob ,n-1)); + sigmaci(2,1)=sigmahat*sqrt((n-1)/chi2inv(1.-prob,n-1)); + catch me + muci(1,1) =muhat; + muci(2,1) =muhat; + sigmaci(1,1)=sigmahat; + sigmaci(2,1)=sigmahat; + end + end + +end + Index: /issm/trunk-jpl/src/m/miscellaneous/norminv_issm.m =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/norminv_issm.m (revision 14182) +++ /issm/trunk-jpl/src/m/miscellaneous/norminv_issm.m (revision 14182) @@ -0,0 +1,9 @@ +% +% wrapper for norminv to avoid using the matlab statistics toolbox. +% +function [x]=norminv_issm(p,mu,sigma) + + x=mu+sigma*sqrt(2.)*erfinv(2.*p-1.); + +end + Index: /issm/trunk-jpl/src/m/qmu/dakota_cdfs.m =================================================================== --- /issm/trunk-jpl/src/m/qmu/dakota_cdfs.m (revision 14181) +++ /issm/trunk-jpl/src/m/qmu/dakota_cdfs.m (revision 14182) @@ -247,5 +247,5 @@ irow=irow+1; cdf(irow,1)=resp(i); - cdf(irow,2)=normcdf(resp(i),mu,sigma); + cdf(irow,2)=normcdf_issm(resp(i),mu,sigma); cdf(irow,3)=(mu-resp(i))/sigma; cdf(irow,4)=(mu-resp(i))/sigma; @@ -256,8 +256,8 @@ for i=1:length(prob) irow=irow+1; - cdf(irow,1)=norminv(prob(i),mu,sigma); + cdf(irow,1)=norminv_issm(prob(i),mu,sigma); cdf(irow,2)=prob(i); - cdf(irow,3)=-norminv(prob(i),0,1); - cdf(irow,4)=-norminv(prob(i),0,1); + cdf(irow,3)=-norminv_issm(prob(i),0,1); + cdf(irow,4)=-norminv_issm(prob(i),0,1); end @@ -267,5 +267,5 @@ irow=irow+1; cdf(irow,1)=mu-sigma*rel(i); - cdf(irow,2)=normcdf(-rel(i),0,1); + cdf(irow,2)=normcdf_issm(-rel(i),0,1); cdf(irow,3)=rel(i); cdf(irow,4)=rel(i); @@ -277,5 +277,5 @@ irow=irow+1; cdf(irow,1)=mu-sigma*grel(i); - cdf(irow,2)=normcdf(-grel(i),0,1); + cdf(irow,2)=normcdf_issm(-grel(i),0,1); cdf(irow,3)=grel(i); cdf(irow,4)=grel(i); Index: /issm/trunk-jpl/src/m/qmu/dakota_moments.m =================================================================== --- /issm/trunk-jpl/src/m/qmu/dakota_moments.m (revision 14181) +++ /issm/trunk-jpl/src/m/qmu/dakota_moments.m (revision 14182) @@ -133,8 +133,15 @@ sigma=std(samp); - muci(1) =mu-tinv(prob,nsamp-1)*sigma/sqrt(nsamp); - muci(2) =mu+tinv(prob,nsamp-1)*sigma/sqrt(nsamp); - sigmaci(1)=sigma*sqrt((nsamp-1)/chi2inv(prob ,nsamp-1)); - sigmaci(2)=sigma*sqrt((nsamp-1)/chi2inv(1-prob,nsamp-1)); + try + muci(1,1) =mu-tinv(prob,nsamp-1)*sigma/sqrt(nsamp); + muci(2,1) =mu+tinv(prob,nsamp-1)*sigma/sqrt(nsamp); + sigmaci(1,1)=sigma*sqrt((nsamp-1)/chi2inv(prob ,nsamp-1)); + sigmaci(2,1)=sigma*sqrt((nsamp-1)/chi2inv(1-prob,nsamp-1)); + catch me + muci(1,1) =mu; + muci(2,1) =mu; + sigmaci(1,1)=sigma; + sigmaci(2,1)=sigma; + end end Index: /issm/trunk-jpl/src/m/qmu/dakota_out_parse.m =================================================================== --- /issm/trunk-jpl/src/m/qmu/dakota_out_parse.m (revision 14181) +++ /issm/trunk-jpl/src/m/qmu/dakota_out_parse.m (revision 14182) @@ -198,5 +198,5 @@ %dstddev=std (data,0); [dmean,dstddev,dmeanci,dstddevci]=... - normfit(data,0.05); + normfit_issm(data,0.05); else dmean =zeros(1,size(data,2)); @@ -206,5 +206,5 @@ for i=1:size(data,2) [dmean(1,i),dstddev(1,i),dmeanci(:,i),dstddevci(:,i)]=... - normfit(data(:,i),0.05); + normfit_issm(data(:,i),0.05); end end Index: /issm/trunk-jpl/src/m/qmu/plot/plot_cdf.m =================================================================== --- /issm/trunk-jpl/src/m/qmu/plot/plot_cdf.m (revision 14181) +++ /issm/trunk-jpl/src/m/qmu/plot/plot_cdf.m (revision 14182) @@ -139,7 +139,7 @@ if strncmpi(cplot,'p',1) && ... exist('ynorm','var') && strncmpi(ynorm,'y',1) - ydata(1:lcdfr(i),i)=norminv(dresp(i).cdf(:,iplot),0,1); + ydata(1:lcdfr(i),i)=norminv_issm(dresp(i).cdf(:,iplot),0,1); else - ydata(1:lcdfr(i),i)= dresp(i).cdf(:,iplot); + ydata(1:lcdfr(i),i)= dresp(i).cdf(:,iplot); end end @@ -237,5 +237,5 @@ end - tick = norminv(yprob,0,1); + tick = norminv_issm(yprob,0,1); set(ax1,'YTick',tick,'YTickLabel',label); ylim([tick(1) tick(end)]) Index: /issm/trunk-jpl/src/m/qmu/plot/plot_hist_norm.m =================================================================== --- /issm/trunk-jpl/src/m/qmu/plot/plot_hist_norm.m (revision 14181) +++ /issm/trunk-jpl/src/m/qmu/plot/plot_hist_norm.m (revision 14182) @@ -224,9 +224,9 @@ if exist('mu','var') && exist('sigma','var') for i=1:ncol - dbelow(end+1)=normcdf(edges( 1),mu(i),sigma(i)); + dbelow(end+1)=normcdf_issm(edges( 1),mu(i),sigma(i)); dhistc(1:size(dhistc,1)-1,end+1)=... - normcdf(edges(2:end ),mu(i),sigma(i))-... - normcdf(edges(1:end-1),mu(i),sigma(i)); - dabove(end+1)=norminv(edges(end),mu(i),sigma(i)); + normcdf_issm(edges(2:end ),mu(i),sigma(i))-... + normcdf_issm(edges(1:end-1),mu(i),sigma(i)); + dabove(end+1)=norminv_issm(edges(end),mu(i),sigma(i)); if exist('descr','var') descr(end+1)={[descr{i} ' norm']}; Index: /issm/trunk-jpl/src/m/qmu/plot/plot_hist_norm_ci.m =================================================================== --- /issm/trunk-jpl/src/m/qmu/plot/plot_hist_norm_ci.m (revision 14181) +++ /issm/trunk-jpl/src/m/qmu/plot/plot_hist_norm_ci.m (revision 14182) @@ -129,5 +129,5 @@ % calculate 95% confidence intervals (same as dakota) [mu(i),sigma(i),muci(:,i),sigmaci(:,i)]=... - normfit(sampr(:,i),0.05); + normfit_issm(sampr(:,i),0.05); end display('Using calculated normal fits.') @@ -191,9 +191,9 @@ if exist('mu','var') && exist('sigma','var') for i=1:ncol - dbelow(end+1)=normcdf(edges( 1),mu(i),sigma(i)); + dbelow(end+1)=normcdf_issm(edges( 1),mu(i),sigma(i)); dhistc(1:size(dhistc,1)-1,end+1)=... - normcdf(edges(2:end ),mu(i),sigma(i))-... - normcdf(edges(1:end-1),mu(i),sigma(i)); - dabove(end+1)=norminv(edges(end),mu(i),sigma(i)); + normcdf_issm(edges(2:end ),mu(i),sigma(i))-... + normcdf_issm(edges(1:end-1),mu(i),sigma(i)); + dabove(end+1)=norminv_issm(edges(end),mu(i),sigma(i)); if exist('descr','var') descr(end+1)={[descr{i} ' norm']}; @@ -204,9 +204,9 @@ if exist('muci','var') && exist('sigmaci','var') for i=1:ncol - dbelow(end+1)=normcdf(edges( 1),muci(1,i),sigmaci(2,i)); + dbelow(end+1)=normcdf_issm(edges( 1),muci(1,i),sigmaci(2,i)); dhistc(1:size(dhistc,1)-1,end+1)=... - normcdf(edges(2:end ),muci(1,i),sigmaci(2,i))-... - normcdf(edges(1:end-1),muci(1,i),sigmaci(2,i)); - dabove(end+1)=norminv(edges(end),muci(1,i),sigmaci(2,i)); + normcdf_issm(edges(2:end ),muci(1,i),sigmaci(2,i))-... + normcdf_issm(edges(1:end-1),muci(1,i),sigmaci(2,i)); + dabove(end+1)=norminv_issm(edges(end),muci(1,i),sigmaci(2,i)); if exist('descr','var') descr(end+1)={[descr{i} ' norm-+']}; @@ -214,9 +214,9 @@ end for i=1:ncol - dbelow(end+1)=normcdf(edges( 1),muci(1,i),sigmaci(1,i)); + dbelow(end+1)=normcdf_issm(edges( 1),muci(1,i),sigmaci(1,i)); dhistc(1:size(dhistc,1)-1,end+1)=... - normcdf(edges(2:end ),muci(1,i),sigmaci(1,i))-... - normcdf(edges(1:end-1),muci(1,i),sigmaci(1,i)); - dabove(end+1)=norminv(edges(end),muci(1,i),sigmaci(1,i)); + normcdf_issm(edges(2:end ),muci(1,i),sigmaci(1,i))-... + normcdf_issm(edges(1:end-1),muci(1,i),sigmaci(1,i)); + dabove(end+1)=norminv_issm(edges(end),muci(1,i),sigmaci(1,i)); if exist('descr','var') descr(end+1)={[descr{i} ' norm--']}; @@ -224,9 +224,9 @@ end for i=1:ncol - dbelow(end+1)=normcdf(edges( 1),muci(2,i),sigmaci(1,i)); + dbelow(end+1)=normcdf_issm(edges( 1),muci(2,i),sigmaci(1,i)); dhistc(1:size(dhistc,1)-1,end+1)=... - normcdf(edges(2:end ),muci(2,i),sigmaci(1,i))-... - normcdf(edges(1:end-1),muci(2,i),sigmaci(1,i)); - dabove(end+1)=norminv(edges(end),muci(2,i),sigmaci(1,i)); + normcdf_issm(edges(2:end ),muci(2,i),sigmaci(1,i))-... + normcdf_issm(edges(1:end-1),muci(2,i),sigmaci(1,i)); + dabove(end+1)=norminv_issm(edges(end),muci(2,i),sigmaci(1,i)); if exist('descr','var') descr(end+1)={[descr{i} ' norm+-']}; @@ -234,9 +234,9 @@ end for i=1:ncol - dbelow(end+1)=normcdf(edges( 1),muci(2,i),sigmaci(2,i)); + dbelow(end+1)=normcdf_issm(edges( 1),muci(2,i),sigmaci(2,i)); dhistc(1:size(dhistc,1)-1,end+1)=... - normcdf(edges(2:end ),muci(2,i),sigmaci(2,i))-... - normcdf(edges(1:end-1),muci(2,i),sigmaci(2,i)); - dabove(end+1)=norminv(edges(end),muci(2,i),sigmaci(2,i)); + normcdf_issm(edges(2:end ),muci(2,i),sigmaci(2,i))-... + normcdf_issm(edges(1:end-1),muci(2,i),sigmaci(2,i)); + dabove(end+1)=norminv_issm(edges(end),muci(2,i),sigmaci(2,i)); if exist('descr','var') descr(end+1)={[descr{i} ' norm++']}; Index: /issm/trunk-jpl/src/m/qmu/plot/plot_normdist_bars.m =================================================================== --- /issm/trunk-jpl/src/m/qmu/plot/plot_normdist_bars.m (revision 14181) +++ /issm/trunk-jpl/src/m/qmu/plot/plot_normdist_bars.m (revision 14182) @@ -126,5 +126,5 @@ if ~isfield(dresp,'mean') || ~isfield(dresp,'stddev') for i=1:length(dresp) - [dresp(i).mean,dresp(i).stddev]=normfit(dresp(i).sample); + [dresp(i).mean,dresp(i).stddev]=normfit_issm(dresp(i).sample); end end @@ -141,5 +141,5 @@ for i=1:length(dresp) descr(i)=cellstr(dresp(i).descriptor); - data(i,:)=norminv(prob,dresp(i).mean,dresp(i).stddev); + data(i,:)=norminv_issm(prob,dresp(i).mean,dresp(i).stddev); end Index: /issm/trunk-jpl/src/m/qmu/plot/plot_rlev_bars_ci.m =================================================================== --- /issm/trunk-jpl/src/m/qmu/plot/plot_rlev_bars_ci.m (revision 14181) +++ /issm/trunk-jpl/src/m/qmu/plot/plot_rlev_bars_ci.m (revision 14182) @@ -133,5 +133,5 @@ [dresp(i).mean,dresp(i).stddev,... dresp(i).meanci,dresp(i).stddevci]=... - normfit(sampr(:,i),0.05); + normfit_issm(sampr(:,i),0.05); display('Using calculated normal fits from sample data.') end @@ -140,7 +140,7 @@ % use minus/plus integer standard deviations sdvect=[-4 -3 -2 -1 0 1 2 3 4]; - dresp(i).cdf(:,2)=normcdf(sdvect,0,1); - dresp(i).cdf(:,1)=norminv(dresp(i).cdf(:,2),... - dresp(i).mean,dresp(i).stddev); + dresp(i).cdf(:,2)=normcdf_issm(sdvect,0,1); + dresp(i).cdf(:,1)=norminv_issm(dresp(i).cdf(:,2),... + dresp(i).mean,dresp(i).stddev); display('Using integer standard deviations for percentages.') @@ -176,5 +176,5 @@ isfield(dresp(i),'stddev') && ~isempty(dresp(i).stddev) descr(end+1)=cellstr([dresp(i).descriptor ' norm']); - cdfr(end+1,:)=norminv(dresp(i).cdf(:,2),dresp(i).mean,dresp(i).stddev); + cdfr(end+1,:)=norminv_issm(dresp(i).cdf(:,2),dresp(i).mean,dresp(i).stddev); end if isfield(dresp(i),'meanci' ) && ~isempty(dresp(i).meanci ) && ... @@ -184,8 +184,8 @@ descr(end+1)=cellstr([dresp(i).descriptor ' norm+-']); descr(end+1)=cellstr([dresp(i).descriptor ' norm++']); - cdfr(end+1,:)=norminv(dresp(i).cdf(:,2),dresp(i).meanci(1),dresp(i).stddevci(2)); - cdfr(end+1,:)=norminv(dresp(i).cdf(:,2),dresp(i).meanci(1),dresp(i).stddevci(1)); - cdfr(end+1,:)=norminv(dresp(i).cdf(:,2),dresp(i).meanci(2),dresp(i).stddevci(1)); - cdfr(end+1,:)=norminv(dresp(i).cdf(:,2),dresp(i).meanci(2),dresp(i).stddevci(2)); + cdfr(end+1,:)=norminv_issm(dresp(i).cdf(:,2),dresp(i).meanci(1),dresp(i).stddevci(2)); + cdfr(end+1,:)=norminv_issm(dresp(i).cdf(:,2),dresp(i).meanci(1),dresp(i).stddevci(1)); + cdfr(end+1,:)=norminv_issm(dresp(i).cdf(:,2),dresp(i).meanci(2),dresp(i).stddevci(1)); + cdfr(end+1,:)=norminv_issm(dresp(i).cdf(:,2),dresp(i).meanci(2),dresp(i).stddevci(2)); end end