How can I evaluate characteristic functions in MatLab?

18 visualizaciones (últimos 30 días)
Martim Zurita
Martim Zurita el 3 de Mzo. de 2021
Editada: Martim Zurita el 15 de Abr. de 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.

Respuesta aceptada

Paul
Paul el 4 de Mzo. de 2021
Editada: Paul el 4 de Mzo. de 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 comentarios
Martim Zurita
Martim Zurita el 4 de Mzo. de 2021
Editada: Martim Zurita el 5 de Mzo. de 2021
Paul, thanks so much for your answer. Not only it solved my problem, it also deepened my knowledge of programming and Fourier analysis.
Paul
Paul el 5 de Mzo. de 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.

Iniciar sesión para comentar.

Más respuestas (1)

Jeff Miller
Jeff Miller el 3 de Mzo. de 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 comentario
Martim Zurita
Martim Zurita el 4 de Mzo. de 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.

Iniciar sesión para comentar.

Productos


Versión

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by