Help understanding why quadgk/quadl does not work

6 次查看（过去 30 天）
-unknown
-unknown 2016-2-12
编辑： -unknown 2016-2-12
Im trying to integrate a specific function, but it is giving me a hard time getting it to work. I want to integrate a product of functions. Im using the build in Regularized Incomplete Beta Function betainc and an own function to calculate the PDF/CDF of a Generalized Hyperbolic Distribution. For completeness of my question, this is the code for the PDF:
function y=ghyppdf(t,lam,alpha,beta,delta,mu)
%Generalized Hyperbolic probability density function (pdf).
% Y=GHYPPDF(X,LAM,ALPHA,BETA,DELTA,MU) returns the pdf of the generalized
% hyperbolic distribution with parameter LAM,
% shape parameter ALPHA, skewness parameter BETA,
% scale parameter DELTA and location parameter MU, evaluated at the
% values in X.
if (delta>0)&(alpha>abs(beta))
kappa = (alpha^2-beta^2)^(lam/2)/(sqrt(2*pi)*alpha^(lam-1/2)*delta^lam*besselk(lam,delta*sqrt(alpha^2-beta^2)));
y = kappa*(delta^2+(t-mu).^2).^((lam-1/2)/2).*besselk(lam-1/2,alpha*sqrt(delta^2+(t-mu).^2)).*exp(beta*(t-mu));
else
y = Inf*ones(size(t));
end;
y(isnan(y))=0;
end
And the CDF:
function y=ghypcdf(x,lam,alpha,beta,delta,mu,starti)
%HYPCDF Hyperbolic cumulative distribution function (cdf).
% Y=HYPCDF(X,ALPHA,BETA,DELTA,MU,STARTI) returns the cdf of the hyperbolic
% distribution with shape parameter ALPHA, skewness parameter BETA,
% scale parameter DELTA and location parameter MU, evaluated at the
% values in X. STARTI is an optional parameter specifying the starting
% point for the integration scheme.
% Convert to a column vector
x = x(:);
% Find the starting point for the integration scheme
feps = 1e-10;
if nargin < 7
starti = mu+beta*delta/sqrt(alpha^2-beta^2)*besselk(lam+1,delta*sqrt(alpha^2-beta^2))/besselk(lam,delta*sqrt(alpha^2-beta^2));
starti = min(starti,min(x))-1;
while ghyppdf(starti,lam,alpha,beta,delta,mu)>feps
starti = starti-1;
end;
end;
n = length(x);
y = zeros(n,1);
[x,ind] = sort(x);
ind = sortrows([ind,(1:n)'],1);
ind = ind(:,2);
x = [starti;x];
warning off MATLAB:quadl:MinStepSize
% Integrate the hyperbolic pdf
for i=1:n
y(i) = quadl('ghyppdf',x(i),x(i+1),[],[],lam,alpha,beta,delta,mu);
end;
warning on MATLAB:quadl:MinStepSize
y = cumsum(y);
y = y(ind);
y(y<0) = 0;%security
y(y>1) = 1;
y(isnan(y)) = 0;%security
end
Now im trying to calculate the following integral in three ways:
1) EV = integral( @(a)betainc(1-ghypcdf(a,lam,alpha,beta,delta,mu),T-Tp,Tp).*a.*ghyppdf(a,lam,alpha,beta,delta,mu),-Inf,Inf,'ArrayValued',true)
2) EV = quadgk( @(a)betainc(1-ghypcdf(a,lam,alpha,beta,delta,mu),T-Tp,Tp).*a.*ghyppdf(a,lam,alpha,beta,delta,mu),-Inf,Inf)
3) EV = quadl( @(a)betainc(1-ghypcdf(a,lam,alpha,beta,delta,mu),T-Tp,Tp).*a.*ghyppdf(a,lam,alpha,beta,delta,mu),-Inf,Inf)
For some reason only the first method works while the other two give me errors saying that the matrix dimensions do not agree. Altough method 1 is working, i want to understand what is going wrong with the others. Based on this i have two questions:
1) Why do i need to include the 'ArrayValued',true part in option 1? It looks like matlab is treating something as an vector/matrix, but it is clearly a scalar function?
2) Why are method 2 and 3 not working, while i explicitly use .* to make a matrix product?
Thanks in advance!

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by