Azzera filtri
Azzera filtri

Can't tell if exprnd is working correctly

11 visualizzazioni (ultimi 30 giorni)
Umberto
Umberto il 1 Lug 2023
Modificato: Umberto il 3 Lug 2023
I'm working on a project involving exponential/Erlang random variables, and I'm trying to validate a model via experiments using exprnd. I'm running into a problem where the results of my code depend in the number of samples I generate using exprnd at one time, no matter how long I run the simulation. I've explored every possibility, and in the end I started wondering whether something might be off with exprnd.
If my understanding is correct, the array exprnd(1/lambda,1,10000) and the array made of 100 copies of exprnd(1/lambda,1,100) should be statistically identical, but when I actually test this via code, there seems to be a difference (but I can't really tell for sure). I'm testing this by simply summing the two vectors. Here's the code I'm using. I'm plotting the difference between the running average of the two sums, which should be zero.
D1=0;
D2=0;
count=0;
lambda=24;
while 1
count=count+1;
d1=0;
for i=1:100
d1=d1+sum(exprnd(1/lambda,1,100));
end
d2=sum(exprnd(1/lambda,1,10000));
D1=D1+d1;
D2=D2+d2;
if rand<.0001
[D1-D2]/count
plot([0 count],[0 0],'r--');hold on;
scatter(count,[D1-D2]/count,'b');
pause(.001);
end
end
Here's a sample plot from the above code.
Right off the bat, it seems strange that the error is consistently negative; I would have thought it should be equally distributed around zero. Also, I can't tell whether or not the error is converging to zero.
Any insight on this would be appreciated.

Risposte (2)

David Goodmanson
David Goodmanson il 1 Lug 2023
Modificato: David Goodmanson il 2 Lug 2023
Hi Umberto,
I believe the issue happens because you are plotting the cumulative variables D1 and D2. Once the difference D1-D2 ends up significantly on one side or the other of zero, it takes awhile to switch back. On one run I had, D1-D2 stayed positive rather than negative for a long time until I got tired of waiting for something to happen.
If you change to plotting noncumulative variables
scatter(count,[d1-d2],'b');
you will see d1-d2 hopping around on both sides of the zero line.
  1 Commento
Umberto
Umberto il 2 Lug 2023
I think you might be right. Also, a t-test fails to reject the null hypothesis of a zero average at any decent confidence level. I think I should keep looking at my code for other errors then.

Accedi per commentare.


Paul
Paul il 2 Lug 2023
Hi Umberto,
I modified the code so that the actual figure of merit that is computed is the difference between the running sample means by including a factor of 1e4 in the denominator.
I hypothesize that the effect you're seeing is just based on the draws. In the code below, for seed = 1, I'm going to guess that exprnd just happened to draw a large outlier (or two or three or ...) on the d1 and that outlier (or outliers?) biases D1 - D2 in favor of D1. That doesn't seem to happen with seed = 3, and maybe goes the other way with seed = 7, though not as much as with seed = 1.
lambda=24;
maxcount = 10000;
figure
hold on
for seed = [1 3 7]
D1=0;
D2=0;
count = 0;
rng(seed)
fom = zeros(maxcount,1);
while count < maxcount
count=count+1;
d1=0;
for i=1:100
d1=d1+sum(exprnd(1/lambda,1,100));
end
d2=sum(exprnd(1/lambda,1,10000));
D1=D1+d1;
D2=D2+d2;
fom(count) = (D1 - D2)/(count*1e4);
%{
if rand<1
[D1-D2]/count/10000;
plot([0 count],[0 0],'r--');hold on;
scatter(count,[D1-D2]/(count*1e4),'b');
%pause(.001);
end
%}
end
plot(fom,'DisplayName',"seed = " + num2str(seed)),grid
end
legend
ylim([-1 1]*1e-4)
Also, those differences are all 2-3 orders of magnitude smaller than expected value of distribution
1/lambda
ans = 0.0417

Prodotti


Release

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by