54 views (last 30 days)

Show older comments

I want to numerically evaluate characteristic functions (CF) of PDFs (definition: https://en.wikipedia.org/wiki/Characteristic_function_(probability_theory) )

For example, the CF of the Normal distribution is

where t is the CF variable, μ is the distribution mean, is its variance and i is the imaginary unit. If I numerically construct a normal distribution as below, how could I numerically estimate its charactertic function?

avg = 1; % average

s = 1; % standard deviation

dx = 0.01;

x = (avg - 5*s):dx:(avg+5*s);

gaussian = 1/sqrt(2*pi*s^2) * exp ( - 0.5 *( x - avg ).^2 / s^2 );

EDIT: I have found an answer that works if one has access to the data generated by the PDF. In the case of a Gaussian distribution with average avg and standard deviation s, one could estimate the CF by the Empirical Characteristic Function (ECF, detailed in Feuerverger & Mureika 1977):

In code,

N = 1e4; % Number of data points

y = m + s*randn(1,N); % N gaussian random numbers with average avg and std s

u = 0.1:0.05:3; % Argument of the characteristic function

ECF = zeros(1,length(u)); % Emperical characteristic function

for j = 1:length(u)

ECF(j) = 1/N * sum( exp ( 1i * u(j) * y ) );

end

The answer provided by Paul below yields equivalent results.

Paul
on 4 Mar 2021

Edited: Paul
on 4 Mar 2021

The CF is the Fourier transform of the pdf. Approximate the pdf as a sum of line segments. Take the sum of the Fourier transforms of each segment.

syms g(x) a b c m w

assume([x a b c m w],'real');

g(x,a,b,c,m) = rectangularPulse(a,b,x)*(m*(x-a)+c); % equation of line segment

G(w,a,b,c,m) = simplify(conj(fourier(g(x,a,b,c,m)))); % Fourier transform of line segment, conj required when using the default FourierParameters

Gfunc = matlabFunction(G)

%{

Gfunc =

function_handle with value:

@(w,a,b,c,m)-1.0./w.^2.*(m.*exp(a.*w.*1i)-m.*exp(b.*w.*1i)-c.*w.*exp(a.*w.*1i).*1i+c.*w.*exp(b.*w.*1i).*1i-a.*m.*w.*exp(b.*w.*1i).*1i+b.*m.*w.*exp(b.*w.*1i).*1i)

%}

% example with standard normal distribution

% distribution parameters

mu = 1;

sigma = 1;

% the pdf

x = -4:.1:4;

gpdf = normpdf(x,mu,sigma);

% generate the CF

w = -4:.01:4;

gcfunc = 0*w;

for ii = 1:numel(x)-1

gcfunc = gcfunc + Gfunc(w,x(ii),x(ii+1),gpdf(ii),(gpdf(ii+1)-gpdf(ii))/(x(ii+1)-x(ii)));

end

% true CF for comparison

cftrue = (exp(1i*mu*w - sigma^2*w.^2/2));

% compare

subplot(211);plot(w,real(gcfunc),w(1:10:end),real(cftrue(1:10:end)),'o'),grid

subplot(212);plot(w,imag(gcfunc),w(1:10:end),imag(cftrue(1:10:end)),'o'),grid

Paul
on 5 Mar 2021

Glad you found it helpful.

I don't know how much the tails of a distribuiton, should they exist, contribute to the CF, or how close the spacing needs to be. Maybe there's a way to check some properties of the computed CF to determine if the answer is reasonable.

Jeff Miller
on 3 Mar 2021

It seems like your question contains its own answer:

mu = 1;

sigma = 1;

t = -2.5:0.01:2.5;

cf = exp(i*mu*t - sigma^2*t.^2/2);

plot(t,cf)

Or maybe I'm missing something--do you want to be able to do this for any pdf, even when you don't have a handy expression for the cf?

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

Start Hunting!