MATLAB Answers

How can I evaluate characteristic functions in MatLab?

54 views (last 30 days)
Martim Zurita
Martim Zurita on 3 Mar 2021
Edited: Martim Zurita on 15 Apr 2021
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.

Accepted Answer

Paul
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
  2 Comments
Paul
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.

Sign in to comment.

More Answers (1)

Jeff Miller
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?
  1 Comment
Martim Zurita
Martim Zurita on 4 Mar 2021
Sorry Jeff, I may not have been clear on my question. Hopefully Paul was still able to grasp what I meant. Nevertheless, thanks for your attention.

Sign in to comment.

Categories

Products


Release

R2018a

Community Treasure Hunt

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

Start Hunting!