MATLAB: Computing the convolution of the difference between random variables

convolutionrandom variables

I have the empirical distributions of Z, and X, where Z= X + Y. I've fit pearson type 4 distribution to X and Z.

What is the best way in matlab to solve for the pdf of Y= Z-X ? I've discretized the pdf of z and x and run a double summation to come up with a vector of y's. However, the result just doesnt seem quite right..

Thanks!

Best Answer

For a first check, I'd run a Monte-Carlo simulation with the empirical CDFs of X and Z.

Best wishes

Torsten.

Related Solutions

MATLAB: How to go about finding the standard normal probability based on the z-score

doc normcdf
doc normpdf

When you know what you want but not sure the name, try something like

>> lookfor normal
realmin                        - Smallest positive normalized floating point number.
randn                          - Normally distributed pseudorandom numbers.
sprandn                        - Sparse normally distributed random matrix.
surfnorm                       - Surface normals.
isonormals                     - Isosurface normals.
cde                            - cd elliptic function with normalized complex argument.
sne                            - sn elliptic function with normalized complex argument.
addfreqcsmenu                  - Add a cs menu to switch between linear and normalized frequency
convertfrequnits               - converts between Normalized, Hz, kHz, etc
histfit                        - Histogram with superimposed fitted normal density.
jbtest                         - Jarque-Bera hypothesis test of composite normality.
lhsnorm                        - Generate a latin hypercube sample with a normal distribution
logncdf                        - Lognormal cumulative distribution function (cdf).
lognfit                        - Parameter estimates and confidence intervals for lognormal data.
logninv                        - Inverse of the lognormal cumulative distribution function (cdf).
lognlike                       - Negative log-likelihood for the lognormal distribution.
lognpdf                        - Lognormal probability density function (pdf).
lognrnd                        - Random arrays from the lognormal distribution.
lognstat                       - Mean and variance for the lognormal distribution.
mvncdf                         - Multivariate normal cumulative distribution function (cdf).
mvnpdf                         - Multivariate normal probability density function (pdf).
mvnrnd                         - Random vectors from the multivariate normal distribution.
normcdf                        - Normal cumulative distribution function (cdf).
normfit                        - Parameter estimates and confidence intervals for normal data.
norminv                        - Inverse of the normal cumulative distribution function (cdf).
normlike                       - Negative log-likelihood for the normal distribution.
normpdf                        - Normal probability density function (pdf).
normplot                       - Displays a normal probability plot.
normrnd                        - Random arrays from the normal distribution.
normspec                       - Plots normal density between specification limits.
normstat                       - Mean and variance for the normal distribution.
logn3fit                       - Fit a 3-param lognormal dist'n using cumulative probabilities.
wgtnormfit                     - Fitting example for a weighted normal distribution.
wgtnormfit2                    - Fitting example for a weighted normal distribution (log(sigma) parameterization).
>>

Judicious search terms help but seeing the list of things related to "normal" lets you find the two functions of interest (plus a lot more depending upon which toolboxes are available, maybe) that might be of use/interest...

MATLAB: Uncertainty calculation for given equation

No. It is not possible to do in any automatic fashion. At least, not with any tools that I know of that are publically available. This is a problem known as statistical tolerancing. There are at least three methods that I can think of to solve the problem:

1. Monte Carlo simulation
2. Truncated Taylor series approximations
3. Taguchi style methods, including what we named modified Taguchi methods

1. Simulate the system using random number generation for the variables.

2. Approximate your function using a first order truncated Taylor series approximation. It is now linear, so the problem reduces to a linear combination of random variables with known distributions. If your random variables are assumed to be Gaussian (Normally distributed), then the solution is trivial.

This scheme can also be extended to a second order Taylor series approximation, although the computations get slightly more involved.

3. Taguchi methods reduce to a numerical integration. Modified Taguchi methods (as described by me and a colleague in a few papers) were shown to be equivalent to a Gauss-Hermite numerical integration. They allow you to compute approximations to the moments of the function. You can then use those estimated moments to infer the distribution of the response from its mean, variance, skewness, and Kurtosis.

You can look online under the phrase "statistical tolerancing modified Taguchi', and you should find some papers that describe the techniques.

https://www.jstor.org/stable/1269802?seq=1#page_scan_tab_contents

Taguchi methods even have a Wiki page. I suppose that were I feeling sufficiently industrious, I should add a page for the modified methods that we developed.

In fact, very long ago I did write some tools that do statistical tolerancing. In fact, it is a credit to MATLAB that my roughly 20 year old code worked fine. Sorry, but I cannot put the code into the public domain. Not really worth the effort anyway, since the solution is so easy without it.

Regardless, if I assume that your variables a,b,c are all uniformly distributed, it really is not that hard to do:

N = 1e8;
a = 5 + (rand(N,1)*2-1)*0.2;
b = 3 + (rand(N,1)*2-1)*0.15;
c = 19 + (rand(N,1)*2-1)*0.01;
y = (a+b)./(c.^2);
mean(y)
ans =
        0.0221605698565819
std(y)
ans =
       0.00040004350708596
skewness(y)
ans =
       0.00018202561148938
kurtosis(y)
ans =
           2.3545148004182

In fact, those estimates of the moments are pretty accurate. Using more accurate methods (fairly high order Gauss-Legendre numerical integration) they can be computed as:

      mean: 0.022160670958632
       std: 0.000400054003691638
  skewness: 0.00012265988869227
  kurtosis: 2.35442773957793

I'd trust the first 8 or 10 of those digits as being correct. So the Monte Carlo estimation was quite reasonable.

hist(y,1000)

See from the histogram that the result does not take on any simple distribution. Of course, the computations above were done using the presumption that your variables were uniformly distributed with the indicated limits. If you intended them to be Normally distributed instead, the results would be slightly different.

Best Answer

Related Solutions

MATLAB: How to go about finding the standard normal probability based on the z-score

MATLAB: Uncertainty calculation for given equation

Related Question