Newsgroup: it.hobby.enigmi
Titolo: (MA) Altro giochino di probabilità

Supponendo di estrarre numeri reali in [0,1] con distr. unif.
e di sommarli fino a superare l'1, a quale valore si arriva mediamente?
Oppure trovare la densità di probabilità.
Ciao
Andrea

 e di sommarli fino a superare l'1, a quale valore si arriva mediamente?
Ho seguito una strada un po' tortuosa, forse c'è un modo più semplice.
Dunque: la probabilità di superare 1 con due soli numeri è
evidentemente:
p(2) = int_0^1 dx1 int_(1-x1)^1 dx2 = 1/2
e il numero che in media si ottiene è m(2), con:
m(2) p(2) = int_0^1 dx1 int_(1-x1)^1 dx2 (x1 + x2) = 2/3.
Analogamente, se 1 viene superato con 3 numeri si ottiene:
p(3) = int_0^1 dx1 int_0^(1-x1) dx2 int_(1-x1-x2)^1 dx3 = 1/3;
m(3) p(3) = 
int_0^1 dx1 int_0^(1-x1) dx2 int_(1-x1-x2)^1 dx3 (x1 + x2 + x3) = 11/24.
Calcolando in questo modo un po' di valori di m(k)p(k) si trova che:
m(k) p(k) = (1/2)(k - 1)(3k + 2)/(k + 1)!
(non ho provato a dimostrarlo ma si dovrebbe poter fare per induzione).
Quindi la media cercata è:
media = sum(k=2, infinito) m(k) p(k) = e/2 = 1.35914...
Maurizio

 Maurizio
Ok, il risultato č giusto!
Ciao
Andrea

 Io avevo seguito una strada un po' diversa, ma sostanzialmente simile
Si riesce in qualche modo a trovare la densità di probabilità?
M.

Maurizio Frigeni <Maurizio...@dont.emailme.please> scritto
nell'articolo <20020605232837160315@[151.100.55.87]>...

Io ho fatto così:
Sia p(x)=prob. di fermarsi in 1..x , con 1<x<2
(la densità di probabilità che cerchiamo è p'(x) )
sia dq(x)=prob. di riuscire a fare con un numero qualunque di tiri
il punteggio x..x+dx , con 0<x<1
Allora per ogni x in [1,2]:
p(x)=somma {prob(di arrivare prima o poi in y..y+dy)*prob(di finire fra 1 e
x con un ultimo tiro)} , con y da 0 a 1
Visto che y<1 le probabilità di saltare in [1,x] con un tiro sono x-1,
se x<1+y, altrimenti finiamo in [1,1+y], con prob. y.
p(x)=S dq(y)*prob(di finire nell' int. [1,min(1+y,x)] con un tiro) ,y=0..1
L'integrale si spezza in due, e le prob. che compaiono equivalgono
all'ampiezza dell'intervallo corrispondente:
p(x)=S dq(y)*prob(di fare 1-y.. 1 con un tiro) ,y=0..x-1
+
S dq(y)*prob(di fare 1-y.. x-y con un tiro) ,y=x-1..1
p(x)=S dq(y)*y ,y=0..x-1
+
S dq(y)*(x-1) ,y=x-1..1
Ora troviamo dq(x):
dq(x)= sum prob(di fare x..x+dx con k tiri) , k=1..inf
dq(x)= dx + x*dx+ x^2/2*dx + x^3/6*dx+...
dq(x)= e^x *dx
infatti, con k tiri:
x    x-x1     x-x1-..-xk-1
Sdx1 Sdx2 ... Sdxk =     x^k/k!
0    0        0
Allora rimaniamo con:
p(x)=S e^y*y*dy ,y=0..(x-1)
+
S e^y*(x-1)*dy ,y=(x-1) ..1
p(x)=(x-2)*e^(x-1) + 1
+
(x-1)*(e -e^(x-1))
p(x)= -e^(x-1) + e*x -e+1
E infine: p'(x)=e - e^(x-1)
La media risulta e/2, e una simulazione mostra che l'andamento è questo,
anche se probabilmente c'è una strada più breve.
Ciao
Andrea

Andrea Artesiani propone

Stavo per rispondere che non mi trovavo d'accordo con i vostri
risultati, quando ho capito di aver interpretato male. Per non
sprecare il lavoro propongo il gioco così come l'avevo capito:
Supponendo di estrarre numeri reali in [0,1] con distr. unif. e di
sommarli fino a superare l'1, quanti lanci servono mediamente?
Oppure trovare la densità di probabilità.
Ciao, Silv:o)

 Ciao, Silv:o)
Questa versione è già stata proposta da me e altri in passato, e anche tu
hai contribuito a risolverla, se non ricordo male.
Comunque per chi non avesse seguito il thread di allora, è un
bell'enigma... magari nella versione più generale in cui si sommano numeri
da 0 a 1 e si vuole superare un valore X reale qualunque.
Ciao
Andrea

 e anche tu hai contribuito a risolverla, se non ricordo male.
ah, sì?
..ho riscoperto l'acqua calda, allora. Poco male: mi sono divertito.
E poi, a posteriori, passare dal risultato del problema che ha
generato questo thread al risultato della variante è cosa immediata,
anche senza aver capito nulla o quasi di come ci si è arrivati.
Ciao, Silv:o)

 sommarli fino a superare l'1, quanti lanci servono mediamente?
Usando il risultato (ottenuto come nel mio post precedente):
p(k) = (k-1)/k!
si trova subito:
<n> = sum(k=2, infinito) k p(k) = e.
Maurizio

 <n> = sum(k=2, infinito) k p(k) = e.
Esatto.
Si poteva sfruttare il tuo lavoro in modo ancora più immediato.
Il risultato medio è e/2. Non c'è motivo di pensare che i lanci necessari a
superare la soglia di 1 siano particolari, per cui in media ogni lancio
contribuisce con 1/2. Ergo, ci vogliono e lanci.
Viceversa, con gli e lanci che normalmente servono a superare l'unità si
ottiene mediamente un risultato di e/2.
Ciao, Silv:o)

 da 0 a 1 e si vuole superare un valore X reale qualunque.
Ho provato a pensarci un po' - mi scuso in anticipo se quanto segue è
già stato scritto in passato.
Chiamo p(k, x) dx la probabilità di ottenere, con k estrazioni di numeri
casuali compresi fra 0 e 1, una somma compresa fra x e x + dx.
Evidentemente si ha:
p(0, x) = delta(x),
p(k, x) = integrale(t, 0, 1) [p(k-1, x-t)],
dove delta(x) è delta di Dirac.
Questo permette di ottenere p(k, x) per ogni valore di k, anche se,
essendo la relazione ricorsiva, non si riesce a trovare direttamente il
risultato come funzione esplicita di k.
Andando però a vedere come è fatta la densità di probabilità p(k, x), si
trova che in ogni intervallo del tipo [n, n+1] (con n intero) essa è un
polinomio di grado k-1. 
Con un po' di pazienza si riesce anche ad indovinare la forma esplicita
di p(k, x); chiamato p(n, k, x) il polinomio che coincide con p(k, x)
nell'intervallo [n, n+1] (k>n), si trova:
p(0, k, x) = x^(k-1)/(k-1)!
p(1, k, x) = ( x^(k-1) - k (x-1)^(k-1) )/(k-1)!
p(2, k, x) = ( x^(k-1) - k (x-1)^(k-1) + (k(k-1)/2) (x-2)^(k-1) )/(k-1)!
e in generale:
p(n, k, x) = somma(i=0, n) [(-1)^i C(k,i) (x-i)^(k-1)]/(k-1)!
La probabilità pt(k, X) di superare la soglia X con k estrazioni è
uguale alla probabilità di ottenere con k-1 estrazioni un numero t in
[X-1, X], moltiplicata per la probabilità (t-X+1) di superare X col tiro
successivo, il tutto sommato sui possibili valori di t:
pt(k, X) = integrale(t, X-1, X) [p(k-1, t)(t - X + 1)].
Ritroviamo così, ad esempio:
pt(k, 1) = (k - 1)/k!
che avevo già trovato per altra strada; inoltre:
pt(k, 2) = (2^(k-1) - k) (k - 2)/k!
ecc.
Adesso è facile trovare il numero medio di estrazioni necessarie a
superare la soglia X:
<k> = somma(k=[X]+1, infinito) [k*pt(k, X)].
Ecco alcuni valori, per X intero:
  X      <k>
--------------------------------------
  1       e
  2       e^2 - e
  3       (2e^3 - 4e^2 + e)/2
  4       (6e^4 - 18e^3 + 12e^2 - e)/6
--------------------------------------
Ora mi piacerebbe sapere innanzitutto se i risultati sono corretti
(ammesso che qualcuno abbia avuto la pazienza di leggere fin qui);
inoltre sarei curioso di sapere se c'è un metodo diretto per trovare i
polinomi p(n, k, x), visto che io li ho ricavati tirando ad indovinare
dopo aver calcolato p(k, x) per vari valori di k e x.
Maurizio

 dopo aver calcolato p(k, x) per vari valori di k e x.
Non so, io ho usato un metodo un po' diverso:
Sia f(x) il numero medio di turni affinché la somma raggiunga x.
E' chiaro che:
f(x)=sum P(fare t)*(1+turni a partire da t per fare x) con t=0..1
f(x)=sum P(fare t)*(1+turni a partire da 0 per fare x-t) con t=0..1
f(x)=S (1+f(x-t))*dt con t=0..1
f(x)=1+S f(x-t)*dt con t=0..1
derivando ambo i membri troviamo:
f'(x)=f(x)-f(x-1)    (1)
Allora per:
x<0 => f(x)=0, partiamo dallo zero e ogni valore <0 è già superato
0<x<1 => f(x-1)=0 => f'(x)=f(x), con evidentemente f(0+)=1 => f(x)=e^x
1<x<2 => f'(x)=f(x)-f(x-1), con f(x-1)=e^(x-1) appena trovata, si risolve
l'eq. diff. e imponendo la continuità in 1, oppure ponendo che soddisfi
l'integrale da cui la (1) deriva, si trova:
f(x)=e^x -(x-1)*e^(x-1)
facendo un altro paio di prove con 2<x<3 e 3<x<4 si intuisce (e si dimostra
facilmente sostituendola nella (1)) la seguente espressione generale:
f(x)=sum (-1)^i*(x-i)^i*e^(x-i)/i!  con i=0..n ; valida per n<=x<=n+1
Tracciandone il grafico si scopre che f(x) è asintotica a 2x
Lo deve essere, visto che per x molto grande i numeri da 0 a 1
che estraggo si possono confondere con la loro media, 0.5, per cui
dovrò estrarre x/0.5=2*x numeri.
E' interessante il fatto che anche f(x)=2x soddisfa sia la (1) che
l'integrale da cui la (1) deriva!
Ciao
Andrea

 da 0 a 1 e si vuole superare un valore X reale qualunque.
Ho trovato un'altra strada possibile rispetto al mio messaggio di ieri,
che sfrutta l'idea già esposta da Silvio: il numero medio di estrazioni
necessarie a superare la soglia X è uguale al doppio della somma media
ottenuta, visto che ogni estrazione contribuisce in media 1/2.
Per trovare la somma media, non è necessario conoscere la densità p(k,
x) della probabilità di ottenere somma x con k estrazioni; è sufficiente
conoscere la densità di probabilità q(x) di ottenere x con un numero
qualsiasi di estrazioni: si tratta della densità già introdotta da
Andrea in un messaggio precedente limitatamente al caso in cui x sia in
[0, 1] e che vale in tal caso e^x.
E' possibile tuttavia trovare q(x) nel caso generale, osservando che
vale l'equazione:
q(x) = integrale(t, x-1, x) q(t)
e che inoltre q(x) = e^x per 0 < x < 1.
Per 1 < x < 2, ad esempio, possiamo scrivere:
q(x) = integrale(t, 1, x) q(t) + integrale(t, x-1, 1) e^t  
cioè:
q(x) = integrale(t, 1, x) q(t) + e - e^(x-1).
Questa si trasforma in un'equazione differenziale:
q'(x) = q(x) - e^(x-1)
con la condizione iniziale:
q(1) = e - 1,
che risolta dà:
q(x) = e^x - x e^(x-1),   per 1 < x < 2.
Analogamente si trova:  
q(x) = e^x - x e^(x-1) + (1/2) x (x-2) e^(x-2),   per 2 < x < 3,
q(x) = e^x - x e^(x-1) + (1/2) x (x-2) e^(x-2) + 
       (1/6) x (x-3)^2 e^(x-3),                     per 3 < x < 4,
cioè in definitiva:
q(x) = somma(i, 0, n) (-1)^i/i! x (x-i)^(i-1) e^(x-i),
per x compreso fra n e n+1 interi.
Conoscendo q(x) possiamo trovare subito la densità di probabilità qt(x)
di ottenere x subito dopo aver superato una soglia prefissata X:
qt(x, X) = integrale(t, x-1, X) q(t)
ed infine il valor medio della somma ottenuta è:
s(X) = integrale(x, X-1, X) [x qt(x)].
In questo modo si trova:
s(1) = e/2    (risultato che già conoscevamo)
s(2) = e(e - 1)/2
s(3) = e(2e^2 - 4e + 1)/4
ecc.
Moltiplicando questi risultati per 2, ritroviamo i valori medi per il
numero di estrazioni già trovati nel mio messaggio di ieri.
Maurizio

 Non so, io ho usato un metodo un po' diverso:
Ho visto solo adesso il messaggio: questo metodo coincide praticamente
con quello che ho postato poco fa.
Maurizio