How can I evaluate characteristic functions in MatLab?

조회 수: 20 (최근 30일)
Martim Zurita
Martim Zurita 2021년 3월 3일
편집: Martim Zurita 2021년 4월 15일
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
Paul 2021년 3월 4일
편집: Paul 2021년 3월 4일
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
Martim Zurita
Martim Zurita 2021년 3월 4일
편집: Martim Zurita 2021년 3월 5일
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 2021년 3월 5일
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.

댓글을 달려면 로그인하십시오.

추가 답변 (1개)

Jeff Miller
Jeff Miller 2021년 3월 3일
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
Martim Zurita
Martim Zurita 2021년 3월 4일
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.

댓글을 달려면 로그인하십시오.

제품


릴리스

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by