dsodf.m

这是一个国外的DST的Matlab工具箱

function erg=dsodf(fhandle,intervalnumber,varargin)
%=========================================================================
% dsodf
%=========================================================================
% y=dsodf(fhandle,intervalnumber,varargin) 
% Discretizes使离散 a given inverse probability逆概率 cdf using the outer discretization离散化 method. 
%
% Input: 
% fhandle: Function handle or function name这是Matlab的概率工具箱的函数名，
% intervalnumber: Number of discretized probability samples离散概率样本数
% varargin: Parameters of fhandle函数带的参数
%
% The cdf must have the form cdf(p,param1,param2,...)产生指定概率分布数据的函数
%
% Output:
% y: Sampled inverse probability cdf with focal elements. 采样逆概率与cdf的焦点元素。
% Size: intnr x length(params(1).ds)... length(params(n).ds)
%
% Usage:
% mu=dsstruct([1,1.4,1])
% sigma=dsstruct([1.4,2.7,0.9;1.5,1.8,0.1])
% y=dsodf('norminv',1000,mu,sigma)
% dscdf(y)
%
%=========================================================================
% Reference : Tonon, F. (2004). "Using random set theory to propagate繁殖,
% epistemic 认识学的uncertainty through a mechanical system." Reliability
% Engineering and System Safety 85(1-3): 169-181.
% Copyright (c) Philipp Limbourg, University of Duisburg-Essen
% www.uni-duisburg-essen.de/informationslogistik/
%=========================================================================
mythres=0.999999/fix(intervalnumber);%fix想0方向元整
erg=dsstruct;
erg.ds=zeros(intervalnumber,3);%intervalnumber行，3列的全为零的矩阵
count=1;
params=[varargin{:}];     %所有的参数在一起的行向量
l=length(params);         %参数的个数
doestr=[];
ind=[];
for i=1:l      %遍历所有的参数变量
    siz=size(params(i).ds,1);   %找出params(i).ds矩阵的第一维的维数，如M行N列矩阵时，siz=M
    if (siz>1)
        doestr=[doestr,siz];    %把所有参数变量的第一维大于1的维数记录下来，成为一个行向量
        ind=[ind,i];            %把>1的维数的位置记下来
    end
end
tempdoe=[];
if (~isempty(doestr))
    tempdoe=fullfact(doestr);     %fullfact全因素试验设计，8行doestr列矩阵
else
    tempdoe=1;
end
doe=ones(size(tempdoe,1),l);      %生成size(tempdoe,1)行l列的全为一的矩阵
doe(:,ind)=tempdoe(:,1:length(ind));    
for i=1:size(doe,1)               %doe矩阵的行数
    prob=1;
    for j=1:l                     %遍历所有的参数
        prob=prob*params(j).ds(doe(i,j),3);
    end
    %  for j=mythres+1-0.9999999:mythres:0.9999999
    temp(1:l)=2;
    doe2=fullfact(temp);
    for k=1:size(doe2,1)
        for k2=1:l
            peval{k2}=params(k2).ds(doe(i,k2),doe2(k,k2));
        end
        t1=mythres+1-0.9999999:mythres:0.9999999;
        t2=t1-mythres;
        borders(:,2*k-1)=feval(fhandle,t1',peval{:});
        borders(:,2*k)=feval(fhandle,t2',peval{:});
    end
    erg.ds(count:count+size(borders,1)-1,1)=min(borders,[],2);
    erg.ds(count:count+size(borders,1)-1,2)=max(borders,[],2);
    erg.ds(count:count+size(borders,1)-1,3)=mythres*prob;
    count=count+size(borders,1);
end
erg=dsnorm(erg);
end