Suppose you are given a bunch of points in the form of a vector x (between "a" and "b") and an accompanying density function p(x)>0. For concreteness, consider
x = [0:0.1:1.0]'; px = 0.3 - (x-0.5).^2;
x = [0:0.1:1.0]'; px = 0.3 - (x-0.5).^2;
This density function is shown in the picture above.
We want to create a N*1 vector "z" of quadrature points that are distributed according to the density p(x). We also want a vector "dz" corresponding to the width of strip belonging to a particular quadrature point. As an additional constraint, we assume that we want z(1) = a, and z(N) = b to match the end-points.
The program attached at the end is able to do this relatively efficiently. In fact, I used the following command to generate the picture above:
[z dz] = GridDensity(x,px,11,1);
GNU Octave Code:
%
% PROGRAM:
% Takes in a PDF or density function, and spits out a bunch of points in
% accordance with the PDF
%
% INPUT:
% x = vector of points. It need *not* be equispaced
% px = vector of same size as x: probability distribution or
% density function. It need not be normalized but has to be positive.
% N = Number of points >= 3. The end points of "x" are included
% necessarily
% Pt = Optional argument. If present, then some plotting.
%
% OUTPUT:
% z = Points distributed according to the density
% hz = width of the "intervals" - useful to apportion domain to points
% if you are doing quadrature with the results, for example.
%
% (*) Sachin Shanbhag, March 5, 2012
% GNU Public License
%
function [z h] = GridDensity(x,px,N,Pt)
npts = 100; % can potentially change
xi = linspace(min(x),max(x),npts)'; % reinterpolate on equi-spaced axis
pint = interp1(x,px,xi,'spline'); % smoothen using splines
ci = cumtrapz(xi,pint);
pint = pint/ci(npts);
ci = ci/ci(npts); % normalize ci
alfa = 1/(N-1); % alfa/2 + (N-1)*alfa + alfa/2
zij = zeros(N,1); % quadrature interval end marker
z = zeros(N,1); % quadrature point
z(1) = min(x);
z(N) = max(x);
%
% ci(Z_j,j+1) = (j - 0.5) * alfa
%
beta = [0.5:1:N-1.5]'*alfa;
zij = [z(1); interp1(ci,xi,beta,'spline'); z(N)];
h = diff(zij);
clear beta;
%
% Quadrature points are not the centroids, but rather the center of masses
% of the quadrature intervals
%
beta = [1:1:N-2]'*alfa;
z(2:N-1) = interp1(ci,xi,beta,'spline');
%
% Some plotting if required
%
if(nargin>3)
subplot(2,1,1)
plot(xi,ci,'b-','LineWidth',2,xi,pint,'r-','LineWidth',2);
axis([min(xi),max(xi)]);
legend('CDF','PDF')
title('Visualization')
xlabel('x');
ylabel('CDF(x)/PDF(x)')
subplot(2,1,2)
plot(z,0.5*ones(size(z)),'ro');
axis([min(xi),max(xi)]);
xlabel('z');
hold on;
for i = 2:N
X = [zij(i); zij(i)];
Y = [0; 1];
plot(X,Y,'b')
endfor
hold off;
endif
endfunction
No comments:
Post a Comment