function Y=NolanStablePDF(X1,nT,alpha,beta,gamma1,delta)
nX=length(X1);
if alpha==1
caucy=0;
thta0=pi/2;
X=(X1-(delta+beta*gamma1*log(gamma1)*2/pi))/gamma1;
else
caucy=-beta*tan(alpha*pi/2);
caucy=floor(caucy*10^8)/10^8;
thta0=atan(beta*tan(alpha*pi/2))/alpha;
X=(X1+caucy-delta)/gamma1;
end
for i1=1:nX
if alpha==1
if beta==0
Y(i1)=1/(pi*(1+X(i1)^2));
else
d=pi/nT;
T=-pi/2:d:pi/2;
for j1=10:nT+1-10
IntV(j1)=Vfun(T(j1),thta0,alpha,beta)* ...
exp(-exp(-pi*X(i1)/(2*beta))*Vfun(T(j1),thta0,alpha,beta));
end
sc=trapz(IntV);
nsc=length(sc);
Y(i1)=sc(nsc)*d*exp(-pi*X(i1)/(2*beta))/(2*abs(beta)*gamma1);
end
else
if X(i1)>caucy
d=(pi/2+thta0)/nT;
T=-thta0:d:pi/2;
for j1=10:nT+1-10
IntV(j1)=Vfun(T(j1),thta0,alpha,beta)* ...
exp(-(X(i1)-caucy)^(alpha/(alpha-1))* ...
Vfun(T(j1),thta0,alpha,beta));
end
sc=trapz(IntV);
nsc=length(sc);
Y(i1)=sc(nsc)*d* ...
alpha*(X(i1)-caucy)^(1/(alpha-1))/(pi*abs(alpha-1)*gamma1);
elseif X(i1)==caucy
dh=gamma(1+1/alpha);
Y(i1)=dh*cos(thta0)/ ...
(pi*(1+caucy^2)^(1/(2*alpha))*gamma1);
else
X2=gamma1*(-X(i1))+delta+caucy;
Y(i1)=NolanStablePDF(X2,nT,alpha,-beta,gamma1,delta);
end
end
end
%[X1',Y']
plot(X1,Y)
评论0