如何计算正态分布随机变量函数的期望？Matlab数据收集优化需求

First, let’s fix your original code to properly estimate the expectation (right now it generates a single sample of e^X instead of an estimate of its mean), then modify it to run multiple trials and analyze the results.

Step 1: Rewrite the Function to Estimate Expectation Correctly #

Your current code produces a single exp(X) value, which is a sample from the distribution of e^X—not an estimate of its expectation. To estimate E[e^X], you need to generate multiple independent samples of X, compute exp(X) for each, then take their mean.

Here’s an updated function that lets you specify how many samples to use per trial:

function [sample_mean] = estimate_exp_expectation(mu, sigma, num_samples)
    % Estimate E[e^X] where X ~ N(mu, sigma^2) using num_samples iid samples
    X = normrnd(mu, sigma, num_samples, 1); % Generate column of samples
    exp_X = exp(X);
    sample_mean = mean(exp_X);
end

Step 2: Run Multiple Trials and Collect Results #

To analyze how your estimates vary, run this function many times and store each result. Preallocate your results array to make the code faster and more memory-efficient.

Example code:

% Define parameters
mu = 0;
sigma = 1;
num_samples_per_trial = 1000; % Samples used to compute one estimate
num_trials = 1000; % Number of times to run the estimation

% Preallocate array for results (avoids slow dynamic growth)
estimates = zeros(num_trials, 1);

% Run all trials
for i = 1:num_trials
    estimates(i) = estimate_exp_expectation(mu, sigma, num_samples_per_trial);
end

Step 3: Plot the Histogram of Estimates #

Plot the distribution of your estimates, and add a vertical line for the true analytical expectation (we know E[e^X] = exp(mu + sigma²/2)—for your parameters, this equals sqrt(e) ≈ 1.6487).

figure;
histogram(estimates, 'Normalization', 'probability'); % Show relative frequencies
hold on;
true_value = exp(mu + sigma^2 / 2);
xline(true_value, 'r--', 'True E[e^X]', 'LineWidth', 2);
xlabel('Estimated Value of E[e^X]');
ylabel('Probability');
title('Distribution of Estimates for E[e^X] (X ~ N(0,1))');
legend('Sample Estimates', 'True Value');
hold off;

Optimization Tips #

• Vectorize for Speed: Skip the trial loop entirely by generating all samples at once. This is far faster in MATLAB for large numbers of trials:

  X = normrnd(mu, sigma, num_samples_per_trial, num_trials); % Matrix of samples
  exp_X = exp(X);
  estimates = mean(exp_X, 1)'; % Column vector of trial means
  

• Preallocate Arrays: Always preallocate arrays (like estimates = zeros(num_trials,1)) instead of letting them grow in loops—this reduces memory overhead and speeds up execution.
• Validate with Analytical Result: Use the closed-form solution exp(mu + sigma²/2) to check that your estimates cluster around the true value.
• Adjust Sample Size: If your histogram is too spread out, increase num_samples_per_trial—more samples per trial will reduce the variance of your estimates.

Bonus: Histogram of Individual e^X Samples #

If you want to visualize the distribution of e^X itself (not the estimates of its mean), generate a large number of samples directly:

large_num_samples = 100000;
X = normrnd(mu, sigma, large_num_samples, 1);
exp_X = exp(X);

figure;
histogram(exp_X, 'Normalization', 'pdf'); % Plot probability density
hold on;
% Overlay the true log-normal PDF (since e^X is log-normal)
x = linspace(min(exp_X), max(exp_X), 1000);
pdf_log_normal = (1./(x * sigma * sqrt(2*pi))) .* exp(-(log(x) - mu).^2/(2*sigma^2));
plot(x, pdf_log_normal, 'r--', 'LineWidth', 2);
xlabel('e^X');
ylabel('Probability Density');
title('Distribution of e^X (X ~ N(0,1))');
legend('Sample Distribution', 'True Log-Normal PDF');
hold off;

内容的提问来源于stack exchange，提问作者liveFreeOrπHard