function [meanx, stddiv, p,q,zlow,zupp] = divmeasure(x,w)
%DIVMEASURE Return deviation measures of x with weights w
%   [meanx, stddiv, p,q,zlow,zupp] =  DIVMEASURE(x,w)
%           x data inputs 
%           w data weights: assume equal weights if w is absent
%           meanx is the mean  
%           stddiv is the standard deviation 
%           p is the forward deviation
%           q is the backward deviation
%           zlow is the lower support, -min(x-meanx)
%           zupp is the upper support, -min(x-meanx)
%           
L = length(x);
if nargin == 1
    w = ones(L,1)/L;
else
    w = w(:)/sum(w);
end;%%%% Estimate PQ

x=x(:);
epsilon = 1e-8;
meanx = w'*x;
z = x-meanx;
stddiv = sqrt(sum(w.*z.^2));
zupp = max(z);
zlow =-min(z);

theta1 = 1e-8;
theta4 = 100/zupp;
theta2 =  theta1 + (theta4-theta1)/3;
theta3 =  theta1 + 2*(theta4-theta1)/3;

% Peform bisection search to find p
while theta4-theta1>epsilon
    tmp2 = sqrt(2*log(sum(exp(theta2*z).*w))/theta2^2);
    tmp3 = sqrt(2*log(sum(exp(theta3*z).*w))/theta3^2);
    %disp([theta2 tmp2 theta3 tmp3])
    if tmp2>tmp3
        theta4 = theta3;
    else
        theta1 = theta2;
    end;
    theta2 =  theta1 + (theta4-theta1)/3;
    theta3 =  theta1 + 2*(theta4-theta1)/3;
end;

p = max((tmp2+tmp3)/2,stddiv);

theta1 = 1e-8;
theta4 = 100/zlow;
theta2 =  theta1 + (theta4-theta1)/3;
theta3 =  theta1 + 2*(theta4-theta1)/3;

% Peform bisection search to find q,
while theta4-theta1>epsilon
    tmp2 = sqrt(2*log(sum(exp(-theta2*z).*w))/theta2^2);
    tmp3 = sqrt(2*log(sum(exp(-theta3*z).*w))/theta3^2);
    %disp([theta2 tmp2 theta3 tmp3])
    if tmp2>tmp3
        theta4 = theta3;
    else
        theta1 = theta2;
    end;
    theta2 =  theta1 + (theta4-theta1)/3;
    theta3 =  theta1 + 2*(theta4-theta1)/3;
end;
q = max((tmp2+tmp3)/2,stddiv);