I looked at simulating Swamis prisoner problem.
I started a review of distributions, exponential and normal being common. How to create a distribution. My Scilab tool has functions to do that, but that is no fun.
Y = k*e^-t/u
The peak value is at t = 0 is k. With u equaling the standard deviation and 5 standard deviations being approximately 99.9% of the values then intuitively k is approximately 5*u.
So I tested it and it worked.
u = 200.
k = 5*u
n = 10000
s = 0
tl = 0
th = k
dt = (th - tl)/n
t = [tl:dt:k]
Numerical integration of the average
for i = 1:n
y(i) = k*%e^(-t(i)/u)
s = s + y(i)
end
a = s/n
Then analytically.
The average value of a function is (1/T)*ʃf(t)dt from t = 0 to T.
Integrating average value a from 0 to T, a = average value.
y = k*e^(t/u)
a (-k*u*e^(-T/u) –-k*u* e^(0/u))/T
a = -k*u*(e^(-T/u) – e^(0/u)/T
With T = k
a = -u*(e^(-T/u) – e^(0/u))
a = -u*(e^(-T/u) – 1)
As T → +inf e^(-T/u) → 0 and for T = k, a = u.
Comparing the Scilab function to what I formulated.
The numerical integration solution and the average value integral solution correlate.
Average via direct integral 99.9955 via numerical integration 100.0455
Scilab comparison.
For using 10 standard deviations in my solution as the peak value, average, median, and cumulative attribution correlate. The standard deviations do not. The problem is a small difference in the distribution between my model and the Scilab function . One curve is a little flatter in the region of the mean value, so the standard deviations are different with Scilab being correct. There is probably a model that is used to ensure the mean and standard deviation are equal. Probably on the net.
Renegades of the mean value, the standard deactivation of my model is always 2x high.
Mine mean 100.0455 std 200.1023
Scilab mean 99.9648 std 99.3087
peaks Scilab 920.8301 Mine 1000.0000
median actual Mine 69.3000 Scilab r 69.3000 calculated 69.3147
Code:
function [cumd,med] = cum_med(y,x)
n = length(y)
cum_sum = 0
y_sum = 0
flag = 0
med = 0
for i = 1:n cum_sum = cum_sum + y(i);end;
for i = 1:n
y_sum = y_sum + y(i)
cumd(i) = 100. * y_sum/cum_sum
if(flag == 0 && cumd(i) >= 50.) then
flag = 1
med = x(i)
end
end
endfunction
function [a,std] = mean_std(y)
n = length(y)
a = 0
for i = 1:n a = a + y(i);end;
a = a /n
s = 0
for i = 1:n s = s + (100. - y(i))^2;end;
std = sqrt(s/n)
endfunction
u = 100.
T = 10. *u // inregration period
k = 10*u // peak value
ai = (-k*u)*(%e^(-T/u) - 1.)/T
n = 10000
dt = T/n
s = 0.
tsum = 0
for i = 1:n
t(i) = tsum
tsum = tsum + dt
y(i) = k*%e^(-t(i)/u)
s = s + y(i)
end
an = s/n
[cd,med] = cum_med(y,t)
mprintf(" Average integral %.4f numerical %.4f\n",ai,an)
[ay,stdy] = mean_std(y)
mprintf(" mean %.4f std %.4f\n",ay,stdy)
r = grand(n,1,"exp",u)
[cdr,medr] = cum_med(y,t)
[asci,stdsci] = mean_std(r)
mprintf(" mean %.4f std %.4f\n",asci,stdsci)
mprintf("peaks %.4f %.4f\n",max(r),max(y))
mprintf("median actual %.4f r %.4f calculated %.4f\n",med,medr,log(2)*u)
r = gsort(r,"g","d")
w1 = scf(1)
clf(w1)
subplot(1,2,1)
plot2d(cd)
plot2d(cdr)
xgrid
subplot(1,2,2)
plot2d(y)
plot2d(r)
xgrid
