How can I evaluate characteristic functions in MatLab?

19 vues (au cours des 30 derniers jours)
Martim Zurita
Martim Zurita le 3 Mar 2021
Modifié(e) : Martim Zurita le 15 Avr 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.

Réponse acceptée

Paul
Paul le 4 Mar 2021
Modifié(e) : Paul le 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 commentaires
Martim Zurita
Martim Zurita le 4 Mar 2021
Modifié(e) : Martim Zurita le 5 Mar 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 le 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.

Connectez-vous pour commenter.

Plus de réponses (1)

Jeff Miller
Jeff Miller le 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 commentaire
Martim Zurita
Martim Zurita le 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.

Connectez-vous pour commenter.

Produits


Version

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by