matlab - 定义代码的X值

我的任务是创建一个脚本，该脚本的行为类似于matlab上的normcdf。

       x=linspace(-5,5,1000); %values for x
       p= 1/sqrt(2*pi) * exp((-x.^2)/2); % THE PDF for the standard normal
       t=cumtrapz(x,p); % the CDF for the standard normal distribution
       plot(x,t); %shows the graph of the CDF

问题是当t值以1:1000而不是-5:5的增量分配时我想知道如何分配正确的x值，即-5:51000到t值输出例如，当我做t(n)时，我得到的结果与normcdf(n)相同。
只是澄清一下：问题是我不能简单地说t(-5)并得到结果=1，就像我在normcdf(1)中所得到的那样，因为cumtrapz计算值被分配给x=1:1000，而不是-5:5。

最佳答案

更新的答案
好的，看了你的评论，下面是你想要做的事情：

x = linspace(-5,5,1000);
p = 1/sqrt(2*pi) * exp((-x.^2)/2);
cdf = cumtrapz(x,p);

q = 3; % Query point
disp(normcdf(q)) % For reference
[~,I] = min(abs(x-q)); % Find closest index
disp(cdf(I)) % Show the value

遗憾的是，没有一个matlab语法可以在一行中很好地完成这项工作，但是如果您将查找最接近的索引抽象到另一个函数中，则可以执行以下操作：

cdf(findClosest(x,q))

function I = findClosest(x,q)
if q>max(x) || q<min(x)
    warning('q outside the range of x');
end
[~,I] = min(abs(x-q));
end

此外，如果您确信查询点q的确切值存在于x，您可以这样做。

cdf(x==q);

但要注意浮点错误你可能会认为某个范围超出了包含某个值的范围，但你几乎不知道它与一个小的舍入erorr不同您可以在下面的示例中看到这一点：

x1 = linspace(0,1,1000); % Range
x2 = asin(sin(x1)); % Ought to be the same thing
plot((x1-x2)/eps); grid on; % But they differ by rougly 1 unit of machine precision

旧答案
据我所知，运行您的代码确实会重现normcdf(x)的结果如果你想像normcdf那样使用erfc。

close all; clear; clc;

x = linspace(-5,5,1000);
cdf = normcdf(x); % Result of normcdf for comparison

%% 1 Trapezoidal integration of normal pd
p = 1/sqrt(2*pi) * exp((-x.^2)/2);
cdf1 = cumtrapz(x,p);

%% 2 But error function IS the integral of the normal pd
cdf2 = (1+erf(x/sqrt(2)))/2;

%% 3 Or, even better, use the error function complement (works better for large negative x)
cdf3 = erfc(-x/sqrt(2))/2;

fprintf('1: Mean error = %.2d\n',mean(abs(cdf1-cdf)));
fprintf('2: Mean error = %.2d\n',mean(abs(cdf2-cdf)));
fprintf('3: Mean error = %.2d\n',mean(abs(cdf3-cdf)));
plot(x,cdf1,x,cdf2,x,cdf3,x,cdf,'k--');

这给了我

1: Mean error = 7.83e-07
2: Mean error = 1.41e-17
3: Mean error = 00 <- Because that is literally what normcdf is doing

如果您的目标不是使用预定义的matlab函数，而是用数值计算结果（即计算错误函数），那么这是一个有趣的挑战，您可以阅读例如here或thisstats stackexchange post作为一个例子，下面的代码通过从第一个链接实现等式2来计算错误函数：

nerf = @(x,n) (-1)^n*2/sqrt(pi)*x.^(2*n+1)./factorial(n)/(2*n+1);

figure(1); hold on;
temp = zeros(size(x)); p =[];
for n = 0:20
    temp = temp + nerf(x/sqrt(2),n);
    if~mod(n,3)
        p(end+1) = plot(x,(1+temp)/2);
    end
end
ylim([-1,2]);
title('\Sigma_{n=0}^{inf}  ( 2/sqrt(pi) ) \times ( (-1)^n x^{2*n+1} ) \div ( n! (2*n+1) )');
p(end+1) = plot(x,cdf,'k--');
legend(p,'n = 0','\Sigma_{n} 0->3','\Sigma_{n} 0->6','\Sigma_{n} 0->9',...
    '\Sigma_{n} 0->12','\Sigma_{n} 0->15','\Sigma_{n} 0->18','normcdf(x)',...
    'location','southeast');
grid on; box on;
xlabel('x'); ylabel('norm. cdf approximations');