function [m,v] = gpstat(k,sigma,theta)
%GPSTAT Mean and variance of the generalized Pareto distribution, optimised for GPU.
%   [M,V] = GPSTAT(K,SIGMA,THETA)
%
%   See also GPSTAT, GPUARRAY.

%   See toolbox/stats/stats/gpstat for implementation details.

%   Copyright 2021 The MathWorks, Inc.

if nargin < 2
    error(message('stats:gpstat:TooFewInputs'));
end
if nargin < 3 || isempty(theta), theta = 0; end

err = internal.stats.statsizechk(3,k,sigma,theta);
if err > 0
    error(message('stats:gpstat:InputSizeMismatch'));
end

outProto = ones("like", internal.stats.dominantType(k,sigma,theta));

% Process each element independently using arrayfun. This causes the
% function to be compiled into an optimised GPU kernel.
[m,v] = arrayfun(@gpstat_scalar, k, sigma, theta, outProto);
end


function [m,v] = gpstat_scalar(k, sigma, theta, outProto)

% Return NaN for out of range parameters.
if sigma < 0
    sigma = NaN("like",sigma);
end

% Calculate mean
if abs(k) < eps
    % k==0 case
    tm = ones("like",k);
elseif k < 1
    tm = 1 ./ (1 - k);
elseif k >= 1
    tm = inf("like", k);
else % k < 0
    tm = nan("like", k);
end

% Calculate variance
if abs(k) < eps
    tv = ones("like",k);
elseif k < 1/2
    tv = 1 ./ ((1-k).^2 .* (1-2.*k));
elseif k >= 1/2
    tv = inf("like",k);
else % k < 0
    tv = nan("like",k);
end

m = theta + sigma .* cast(tm,"like",outProto);
v = sigma.^2 .* cast(tv,"like",outProto);
end