function q=quantile(p,x,w);
% outputs three estimates of the quantile xp, 0<p<1. x are the sample values
% and w are the weights attached to them. Assume weights (approx) add to
% one.  The first method ignores the weights (assumes all equal).
% linearly nterpolates between (xn,F(xn))  and (x_n+1 F(x_n+1))
v=ones(1,length(x)); v=v/sum(v);
[x,i]=sort(x);
v=v(i);
cdf=cumsum(v);
i=sum(cdf<=p);  % index of largest point such that F(x)<p. 
z=[x(i) x(i+1)]; F=[cdf(i) cdf(i+1)];
wt=(F(2)-p)/(F(2)-F(1));
q=x(i)*wt+x(i+1)*(1-wt);
cdf=cumsum(w/sum(w));
i=sum(cdf<=p);  % index of largest point such that F(x)<p. 
z=[x(i) x(i+1)]; F=[cdf(i) cdf(i+1)];
wt=(F(2)-p)/(F(2)-F(1));
q=[q x(i)*wt+x(i+1)*(1-wt)];