35

Ostatnio przyglądałem się symulacji Monte Carlo i używałem jej do przybliżania stałych, takich jak (okrąg wewnątrz prostokąta, obszar proporcjonalny). $\pi$

Nie jestem jednak w stanie wymyślić odpowiedniej metody aproksymacji wartości $e$ [liczby Eulera] przy użyciu integracji Monte Carlo.

Czy masz jakieś wskazówki, jak to zrobić?

— statystyki newwbie12345
źródło

7

Jest na to wiele sposobów. To, że tak jest, może stać się oczywiste, rozważając, co robi Rpolecenie 2 + mean(exp(-lgamma(ceiling(1/runif(1e5))-1))). (Jeśli przeszkadza Ci korzystanie z funkcji log Gamma, zastąp ją 2 + mean(1/factorial(ceiling(1/runif(1e5))-2)), która używa tylko dodawania, mnożenia, dzielenia i obcinania, i ignoruj ostrzeżenia o przepełnieniu). Bardziej interesujące mogą być wydajne symulacje: czy możesz zminimalizować liczbę kroki obliczeniowe potrzebne do oszacowania

e

$e$ do dowolnej dokładności?

— whuber

4

Cóż za wspaniałe pytanie! Z niecierpliwością czekam na odpowiedzi innych. Jednym ze sposobów, w jaki naprawdę możesz zwrócić uwagę na to pytanie - być może kolejne pół tuzina odpowiedzi - byłoby zrewidowanie pytania i poproszenie o skuteczne odpowiedzi, jak sugeruje whuber. To jest jak kocimiętka dla użytkowników CV.

— Sycorax mówi Przywróć Monikę

1

@EngrStudent Nie jestem pewien, czy geometryczny analog istnieje dla

e

$e$ . Po prostu nie jest to naturalnie (zamierzona gra słów) wielkość geometryczna, taka jak

π

$\pi$ .

— Aksakal

6

@Aksakal

e

$e$ jest wyjątkowo geometryczną wielkością. Na najbardziej elementarnym poziomie pojawia się naturalnie w wyrażeniach dotyczących obszarów związanych z hiperbolami. Na nieco bardziej zaawansowanym poziomie jest ściśle związany ze wszystkimi funkcjami okresowymi, w tym funkcjami trygonometrycznymi, których zawartość geometryczna jest oczywista. Prawdziwym wyzwaniem jest to, że tak łatwo symulować wartości związane z

e

$e$ !

— whuber

2

@StatsStudent:

e

$e$ sam w sobie nie jest interesujący. Jeśli jednak prowadzi to do obiektywnych estymatorów wielkości, takich jak

\exp {\int_{0}^{x} f (y) d G (y)}

$\exp\left\{\int_0^x f(y) \text{d}G(y)\right\}$ może się to okazać najbardziej przydatne w algorytmach MCMC.

— Xi'an

34

W tym artykule opisano prosty i elegancki sposób oszacowania przez Monte Carlo . Artykuł dotyczy nauczania . Dlatego podejście wydaje się idealnie pasować do twojego celu. Pomysł opiera się na ćwiczeniu z popularnego rosyjskiego podręcznika teorii prawdopodobieństwa autorstwa Gnedenko. Patrz przykład 22 na str. 183 $e$ $e$

Zdarza się tak, że , gdzie jest zmienną losową, która jest zdefiniowana w następujący sposób. Jest to minimalna liczba taka, że i są liczbami losowymi z rozkładu równomiernego na . Piękne, prawda ?! $E[\xi]=e$ $\xi$ $n$ $\sum_{i=1}^n r_i>1$ $r_i$ $[0,1]$

Ponieważ jest to ćwiczenie, nie jestem pewien, czy fajnie jest dla mnie opublikować rozwiązanie (dowód) tutaj :) Jeśli chcesz sam to udowodnić, oto wskazówka: rozdział nazywa się „Chwile”, które powinny wskazywać jesteś we właściwym kierunku.

Jeśli chcesz wdrożyć go samodzielnie, nie czytaj dalej!

Jest to prosty algorytm do symulacji Monte Carlo. Narysuj jednolity losowy, a następnie kolejny i tak dalej, aż suma przekroczy 1. Liczba losowanych losów jest twoją pierwszą próbą. Powiedzmy, że masz:

 0.0180
 0.4596
 0.7920

Potem twoja pierwsza próba sprawiła, że 3. Kontynuuj te próby, a zauważysz, że średnio dostajesz . $e$

Poniżej znajduje się kod MATLAB, wynik symulacji i histogram.

N = 10000000;
n = N;
s = 0;
i = 0;
maxl = 0;
f = 0;
while n > 0
    s = s + rand;
    i = i + 1;
    if s > 1
        if i > maxl
            f(i) = 1;
            maxl = i;
        else
            f(i) = f(i) + 1;
        end
        i = 0;
        s = 0;
        n = n - 1;
    end
end

disp ((1:maxl)*f'/sum(f))
bar(f/sum(f))
grid on

f/sum(f)

Wynik i histogram:

2.7183


ans =

  Columns 1 through 8

         0    0.5000    0.3332    0.1250    0.0334    0.0070    0.0012    0.0002

  Columns 9 through 11

    0.0000    0.0000    0.0000

AKTUALIZACJA: Zaktualizowałem swój kod, aby pozbyć się szeregu wyników próbnych, aby nie zajmował pamięci RAM. Wydrukowałem również oszacowanie PMF.

Aktualizacja 2: Oto moje rozwiązanie Excel. Umieść przycisk w programie Excel i połącz go z następującym makrem VBA:

Private Sub CommandButton1_Click()
n = Cells(1, 4).Value
Range("A:B").Value = ""
n = n
s = 0
i = 0
maxl = 0
Cells(1, 2).Value = "Frequency"
Cells(1, 1).Value = "n"
Cells(1, 3).Value = "# of trials"
Cells(2, 3).Value = "simulated e"
While n > 0
    s = s + Rnd()
    i = i + 1
    If s > 1 Then
        If i > maxl Then
            Cells(i, 1).Value = i
            Cells(i, 2).Value = 1
            maxl = i
        Else
            Cells(i, 1).Value = i
            Cells(i, 2).Value = Cells(i, 2).Value + 1
        End If
        i = 0
        s = 0
        n = n - 1
    End If
Wend


s = 0
For i = 2 To maxl
    s = s + Cells(i, 1) * Cells(i, 2)
Next


Cells(2, 4).Value = s / Cells(1, 4).Value

Rem bar (f / Sum(f))
Rem grid on

Rem f/sum(f)

End Sub

Wprowadź liczbę prób, na przykład 1000, w komórce D1 i kliknij przycisk. Oto jak powinien wyglądać ekran po pierwszym uruchomieniu:

AKTUALIZACJA 3: Silverfish zainspirował mnie do innej drogi, nie tak eleganckiej jak pierwsza, ale wciąż fajnej. Obliczył objętości n-simpleksów przy użyciu sekwencji Sobola .

s = 2;
for i=2:10
    p=sobolset(i);
    N = 10000;
    X=net(p,N)';
    s = s + (sum(sum(X)<1)/N);
end
disp(s)

2.712800000000001

Przypadkowo napisał pierwszą książkę o metodzie Monte Carlo, którą przeczytałem w szkole średniej. Moim zdaniem jest to najlepsze wprowadzenie do metody.

AKTUALIZACJA 4:

Silverfish w komentarzach sugerował prostą implementację formuły Excel. Taki wynik uzyskuje się po jego podejściu po około 1 milionie losowych liczb i 185 000 prób:

Oczywiście jest to znacznie wolniejsze niż wdrożenie Excel VBA. Zwłaszcza jeśli zmodyfikujesz mój kod VBA, aby nie aktualizować wartości komórek w pętli i zrobisz to dopiero po zebraniu wszystkich statystyk.

AKTUALIZACJA 5

Xi'an za rozwiązanie nr 3 jest ściśle powiązany (lub nawet taki sam w pewnym sensie jako komentarz na JWG w wątku). Trudno powiedzieć, kto wpadł na ten pomysł jako pierwszy Forsythe lub Gnedenko. Oryginalna edycja Gnedenko z 1950 roku w języku rosyjskim nie zawiera sekcji Problemy w rozdziałach. Tak więc nie mogłem znaleźć tego problemu na pierwszy rzut oka, gdzie jest w późniejszych wydaniach. Może został dodany później lub zakopany w tekście.

Jak skomentowałem w odpowiedzi Xi'ana, podejście Forsythe'a wiąże się z innym interesującym obszarem: rozkładem odległości między pikami (ekstrema) w losowych sekwencjach (IID). Średnia odległość zdarza się wynosić 3. Sekwencja dolna w podejściu Forsythe'a kończy się na dole, więc jeśli będziesz kontynuować próbkowanie, w pewnym momencie dostaniesz kolejne dno, a potem inne. Możesz śledzić odległość między nimi i zbudować rozkład.

— Aksakal
źródło

Wow, to super! Czy możesz dodać akapit lub dwa wyjaśniające, dlaczego to działa?

— Sycorax mówi Przywróć Monikę

7

(+1) Genialne! Odpowiedź zasługuje na najwyższą ocenę, ponieważ opiera się tylko na jednolitych symulacjach. I nie używa żadnego przybliżenia, tylko z powodu Monte Carlo. To, że łączy się ponownie z Gnedenko, jest dodatkowym atutem.

— Xi'an

2

Fajne! Oto kod Mathematica dla tego samego, jako

\dots

$\dots$ Mean[Table[ Length[NestWhileList[(Random[]+#) &, Random[], #<1&]], {10^6}]]

— jednowierszowego

4

@wolfies Następujące bezpośrednie tłumaczenie Rrozwiązania, które zamieściłem w odpowiedzi n=10^6; 1. / Mean[UnitStep[Differences[Sort[RandomReal[{0, n}, n + 1]]] - 1]]

— Xi'ana,

1

Opublikowałem „dlaczego to znaczy ?” pytanie samo w sobie ; Podejrzewam, że moje rozwiązanie szkicu (które natychmiast przyszło mi do głowy jako „oczywista” wizualizacja problemu) niekoniecznie jest tak, jak zamierzali to zrobić rosyjscy studenci! Dlatego alternatywne rozwiązania byłyby bardzo mile widziane.

e

$e$

— Silverfish,

19

Sugeruję, by głosować na odpowiedź Aksakala. Jest bezstronny i opiera się wyłącznie na metodzie generowania odchyleń jednostkowych jednostek.

Moja odpowiedź może być dowolnie sprecyzowana, ale nadal jest stronnicza od prawdziwej wartości . $e$

Odpowiedź Xi'ana jest poprawna, ale myślę, że jej zależność od funkcji lub sposobu generowania losowych odchyleń Poissona jest nieco okrągła, gdy celem jest przybliżenie . $\log$ $e$

Oszacowanie przez ładowanie początkowe $e$

Zamiast tego rozważ procedurę ładowania. Jeden ma dużą liczbę obiektów które są rysowane z zamianą na wielkość próbki . Przy każdym losowaniu prawdopodobieństwo nie rysowania określonego obiektu wynosi i jest takich losowań. Prawdopodobieństwo pominięcia określonego obiektu we wszystkich losowaniach wynosi $n$ $n$ $i$ $1-n^{-1}$ $n$ $p=(1-\frac{1}{n})^n.$

Ponieważ zakładam, że wiemy, że

\exp (- 1) = lim_{n \to \infty} {(1 - \frac{1}{n})}^{n}

$\exp(-1)=\lim_{n\to\infty}\left(1-\frac{1}{n}\right)^n$

więc możemy też napisać

\exp (- 1) \approx \hat{p} = \sum_{i = 1}^{m} \frac{I_{i \in B_{j}}}{m}

$\exp(-1)\approx \hat{p}=\sum_{i=1}^m\frac{\mathbb{I}_{i\in B_j}}{m}$

Oznacza to, że nasze oszacowanie można znaleźć poprzez oszacowanie prawdopodobieństwa pominięcia konkretnej obserwacji w bootstrapie replikuje w wielu takich replikacjach - tj. występowania obiektu w bootstrapie. $p$ $m$ $B_j$ $i$

Są dwa źródła błędów w tym przybliżeniu. Skończone zawsze oznacza, że wyniki są przybliżone, tj. Oszacowanie jest stronnicze. Dodatkowo będzie oscylować wokół prawdziwej wartości, ponieważ jest to symulacja. $n$ $\hat{p}$

Uważam to podejście za urocze, ponieważ student lub inna osoba z wystarczającą ilością rzeczy do zrobienia może przybliżyć za pomocą talii kart, stosu małych kamieni lub innych dostępnych przedmiotów, w tym samym stylu, co osoba może oszacować używając kompasu, prostej krawędzi i drobinek piasku. Myślę, że to fajne, kiedy matematykę można oddzielić od nowoczesnych udogodnień, takich jak komputery. $e$ $\pi$

Wyniki

Przeprowadziłem kilka symulacji dla różnej liczby replik ładowania początkowego. Standardowe błędy są szacowane w normalnych odstępach czasu.

Zauważ, że wybór liczby obiektów ładowanych początkowo ustanawia absolutną górną granicę dokładności wyników, ponieważ procedura Monte Carlo szacuje a zależy tylko od . Ustawienie niepotrzebnie dużej wartości tylko obciąży komputer, albo dlatego, że potrzebujesz tylko „przybliżonego” przybliżenia do albo dlatego, że odchylenie zostanie zatłoczone przez wariancję z powodu Monte Carlo. Te wyniki są dla a jest dokładne z dokładnością do trzeciego miejsca po przecinku. $n$ $p$ $p$ $n$ $n$ $e$ $n=10^3$ $p^{-1}\approx e$

Ten wykres pokazuje, że wybór ma bezpośrednie i głębokie konsekwencje dla stabilności w . Niebieska linia przerywana pokazuje a czerwona linia pokazuje . Zgodnie z oczekiwaniami, zwiększenie wielkości próby daje coraz dokładniejsze szacunki . $m$ $\hat{p}$ $p$ $e$ $\hat{p}$

Napisałem do tego żenująco długi skrypt R. Sugestie dotyczące ulepszeń można przesłać na odwrocie rachunku za 20 USD.

library(boot)
library(plotrix)
n <- 1e3

## if p_hat is estimated with 0 variance (in the limit of infinite bootstraps), then the best estimate we can come up with is biased by exactly this much:
approx <- 1/((1-1/n)^n)

dat <- c("A", rep("B", n-1))
indicator <- function(x, ndx)   xor("A"%in%x[ndx], TRUE) ## Because we want to count when "A" is *not* in the bootstrap sample

p_hat <- function(dat, m=1e3){
    foo <- boot(data=dat, statistic=indicator, R=m) 
    1/mean(foo$t)
} 

reps <- replicate(100, p_hat(dat))

boxplot(reps)
abline(h=exp(1),col="red")

p_mean <- NULL
p_var <- NULL
for(i in 1:10){
    reps <- replicate(2^i, p_hat(dat))
    p_mean[i] <- mean(reps)
    p_var[i] <- sd(reps)
}
plotCI(2^(1:10), p_mean, uiw=qnorm(0.975)*p_var/sqrt(2^(1:10)),xlab="m", log="x", ylab=expression(hat(p)), main=expression(paste("Monte Carlo Estimates of ", tilde(e))))
abline(h=approx, col='red')

— Sycorax mówi Przywróć Monikę
źródło

4

+1 Ma to sens. Czy jest szansa, że możesz podzielić się swoim kodem, jeśli go napisałeś?

— Antoni Parellada,

2

Choć może to być dowolnie dokładne, ostatecznie jest niezadowalająca, ponieważ tylko symuluje zbliżania się zamiast sama.

e

$e$

e

$e$

— whuber

1

Pewnie. Zakończyłbyś się tylko jednym powtórzeniem wewnątrz drugiego, co jest w zasadzie takie samo jak teraz.

— Sycorax mówi Przywróć Monikę

1

@ whuber Naprawdę nie widzę różnicy między arbitralnie dokładnym przybliżeniem do arbitralnie dokładnego przybliżenia do , a arbitralnie dokładnym przybliżeniem do samego .

e

$e$

e

$e$

— jwg

1

@jwg Oprócz tego, że jest to ważne z punktu widzenia koncepcyjnego, jest także praktycznie ważne, ponieważ implementacja aproksymacji do aproksymacji wymaga śledzenia dokładności każdego z tych dwóch aproksymacji. Ale musiałbym się zgodzić, że gdy oba przybliżenia są akceptowalnie dobre, ogólne podejście rzeczywiście jest w porządku.

— whuber

14

Rozwiązanie 1:

Dla Poissona dystrybucja Dlatego, jeśli , co oznacza, że możesz oszacować za pomocą symulacji Poissona. Symulacje Poissona można uzyskać z generatora rozkładu wykładniczego (jeśli nie w najbardziej efektywny sposób). $\mathcal{P}(\lambda)$

P. (X = k) = \frac{λ^{k}}{k!} {mi}^{- λ}

$\mathbb{P}(X=k)=\frac{\lambda^k}{k!}\,e^{-\lambda}$

X \sim P (1)

$X\sim\mathcal{P}(1)$

P. (X = 0) = P. (X = 1) = {mi}^{- 1}

$\mathbb{P}(X=0)=\mathbb{P}(X=1)=e^{-1}$

e^{- 1}

$e^{-1}$

Uwaga 1: Jak omówiono w komentarzach, jest to dość skomplikowany argument, ponieważ symulowanie z rozkładu Poissona lub równoważnie rozkładu wykładniczego może być trudne do wyobrażenia bez udziału log lub funkcji exp ... Ale potem W. Huber doszedł do uratowanie tej odpowiedzi za pomocą najbardziej eleganckiego rozwiązania opartego na zamówionych mundurach. Co jest jednak przybliżeniem , ponieważ rozkład jednolitych odstępów jest Beta , co sugeruje, że który jest zbieżny do jako $U_{(i:n)}-U_{(i-1:n)}$ $\mathfrak{B}(1,n)$
$P. (n {U_{(ja : n)} - U_{(ja - 1 : n)}} \geq 1) = {(1 - \frac{1}{n})}^{n}$ $\mathbb{P}(n\{U_{(i:n)}-U_{(i-1:n)}\}\ge 1)=\left(1-\frac{1}{n}\right)^n$ $e^{-1}$ $n$ rośnie do nieskończoności. Oprócz tego, że odpowiada na komentarze, generator wykładniczy von Neumanna z 1951 r. Używa tylko generacji jednorodnych.

Rozwiązanie 2:

Innym sposobem uzyskania reprezentacji stałej jako całki jest przypomnienie, że gdy to który jest również rozkładem . Dlatego Drugie podejście do przybliżenia przez Monte Carlo ma zatem symulować pary normalne i monitorować częstotliwość razy . W pewnym sensie jest to przeciwieństwo aproksymacji Monte Carlo związane z częstotliwością razy ... $e$

X_{1}, X_{2} \overset{iid}{\sim} N (0, 1)

$X_1,X_2\stackrel{\text{iid}}{\sim}\mathfrak{N}(0,1)$

(X_{1}^{2} + X_{2}^{2}) \sim χ_{1}^{2}

$(X_1^2+X_2^2)\sim\chi^2_1$

E (1 / 2)

$\mathcal{E}(1/2)$

P (X_{1}^{2} + X_{2}^{2} \geq 2) = 1 - {1 - \exp (- 2 / 2)} = e^{- 1}

$\mathbb{P}(X_1^2+X_2^2\ge 2)=1-\{1-\exp(-2/2)\}=e^{-1}$

e

$e$

(X_{1}, X_{2})

$(X_1,X_2)$

X_{1}^{2} + X_{2}^{2} \geq 2

$X_1^2+X_2^2\ge 2$

π

$\pi$

X_{1}^{2} + X_{2}^{2} < 1

$X_1^2+X_2^2<1$

Rozwiązanie 3:

Mój kolega z Uniwersytetu w Warwick, M. Pollock, wskazał inne przybliżenie Monte Carlo zwane metodą Forsythe'a : chodzi o to, by uruchomić sekwencję ujednoliconych generacji aż do . Oczekiwanie na odpowiednią regułę zatrzymania, , która jest liczbą przypadków, w których jednorodna sekwencja spadła, wynosi wtedy podczas gdy prawdopodobieństwo, że jest nieparzyste, wynosi ! ( Metoda Forsythe'a faktycznie ma na celu symulację z dowolnej gęstości formy , stąd jest bardziej ogólna niż przybliżenie i .) $u_1,u_2,...$ $u_{n+1}>u_{n}$ $N$ $e$ $N$ $e^{-1}$ $\exp G(x)$ $e$ $e^{-1}$

Jest to dość równoległe do podejścia Gnedenko zastosowanego w odpowiedzi Aksakala , więc zastanawiam się, czy jedno można wyprowadzić z drugiego. Przynajmniej oba mają ten sam rozkład z masą prawdopodobieństwadla wartości . $1/n!$ $n$

Szybka implementacja metody Forsythego R polega na rezygnacji z dokładnego wykonywania sekwencji mundurów na rzecz większych bloków, co pozwala na równoległe przetwarzanie:

use=runif(n)
band=max(diff((1:(n-1))[diff(use)>0]))+1
bends=apply(apply((apply(matrix(use[1:((n%/%band)*band)],nrow=band),
2,diff)<0),2,cumprod),2,sum)

— Xi'an
źródło

12

Dopóki ktoś wie, jak przeprowadzić symulację Poissona, nie znając .

e

$e$

— Glen_b

5

Jeśli wywołam generator R rpoiss (), mogę udawać, że nie znam . Mówiąc poważniej, można wygenerować zmienne wykładnicze [używając funkcji zamiast ], dopóki suma nie przekroczy a wynikowa liczba minus jeden to Poisson .

e

$e$

E (1)

$\mathcal{E}(1)$

\log

$\log$

e

$e$

1

$1$

P (1)

$\mathcal{P}(1)$

— Xi'an

5

\log

$\log$

\exp

$\exp$ n <- 1e5; 1/mean(n*diff(sort(runif(n+1))) > 1)

3

x_{n}

$x_n$

\sum^{n} x_{i}

$\sum^n x_i$

x_{n}

$x_n$

1 - \sum^{n - 1} x_{i}

$1 - \sum^{n-1} x_i$

1 - \sum^{n} x_{i}

$1 - \sum^{n} x_i$

1 - \sum^{n - 1} x_{i}

$1 - \sum^{n-1} x_i$

3

n + 1

$n+1$

n

$n$

7

To nie rozwiązanie ... tylko szybki komentarz, który jest za długi na pole komentarza.

Aksakal

Aksakal opublikował rozwiązanie, w którym obliczamy oczekiwaną liczbę standardowych rysunków Uniform, które należy pobrać, tak aby ich suma przekroczyła 1. W Mathematica moim pierwszym sformułowaniem było:

mrM := NestWhileList[(Random[] + #) &, Random[], #<1 &]

Mean[Table[Length[mrM], {10^6}]]

EDYCJA: Właśnie się z tym pobawiłem, a następujący kod (ta sama metoda - także w Mma - tylko inny kod) jest około 10 razy szybszy:

Mean[Table[Module[{u=Random[], t=1},  While[u<1, u=Random[]+u; t++]; t] , {10^6}]]

Xian / Whuber

Whuber zasugerował szybki fajny kod do symulacji rozwiązania Xian 1:

Wersja R: n <- 1e5; 1/mean(n*diff(sort(runif(n+1))) > 1)

Wersja Mma: n=10^6; 1. / Mean[UnitStep[Differences[Sort[RandomReal[{0, n}, n + 1]]] - 1]]

który zauważa, że jest 20 razy szybszy niż pierwszy kod (lub około dwa razy szybciej niż nowy kod powyżej).

Dla zabawy pomyślałem, że byłoby interesujące sprawdzić, czy oba podejścia są tak skuteczne (w sensie statystycznym). W tym celu wygenerowałem 2000 oszacowań e przy użyciu:

Metoda Aksakala: dane A.
Metoda Xiana 1 przy użyciu kodu Whuber: dataB

... oba w Mathematica . Poniższy diagram porównuje nieparametryczne oszacowanie gęstości jądra wynikowych zestawów danych dataA i dataB.

Tak więc, chociaż kod Whubera (czerwona krzywa) jest około dwa razy szybszy, metoda nie wydaje się być „niezawodna”.

— wilków
źródło

Pionowa linia w miejscu prawdziwej wartości znacznie poprawiłaby ten obraz.

— Sycorax mówi Przywróć Monikę

1

To bardzo interesujące spostrzeżenie, dziękuję. Ponieważ połowa szerokości będzie skalować się kwadratowo wraz z rozmiarem symulacji, a połowa szerokości metody Xi'ana jest około dwa razy większa niż metody Aksakala, a następnie wykonanie czterech razy większej liczby iteracji sprawi, że będą one równie dokładne. Pozostaje pytanie, ile wysiłku potrzeba w każdej iteracji: jeśli jedna iteracja metody Xi'ana zajmuje mniej niż jedną czwartą wysiłku, wówczas metoda ta byłaby jeszcze bardziej wydajna.

— whuber

1

n

$n$

1

running four times as many iterations will make them equally accurate

10^{6}

$10^6$

10^{6}

$10^6$

1

Dobra robota z kodem - trudno będzie wiele poprawić.

— whuber

2

Metoda wymagająca bezbożnej ilości próbek

$f(x) = e^x$ $\frac{ \bar x \sqrt{12}}{ \sqrt{n}} \dot \sim N(0,1)$ $e$ $\infty$ $\infty$

$N(0,1)$ $e^{x}$

$N(0,1)$ $\hat \phi(x)$ $\phi(\sqrt(2) ) = (2 \pi)^{-1/2} e^{-1}$ $e = \hat \phi(\sqrt{2}) \sqrt{2 \pi}$

$\sqrt{2}$ $\sqrt{2\pi}$

Metoda wymagająca bardzo niewielu próbek, ale powodująca bezbożną ilość błędów numerycznych

Całkowicie głupia, ale bardzo skuteczna odpowiedź na podstawie mojego komentarza:

$X \sim \text{uniform}(-1, 1)$ $Y_n = | (\bar x)^n|$ $\hat e = (1 - Y_n)^{-1/Y_n}$

To zbiegnie się bardzo szybko, ale również napotka ekstremalny błąd numeryczny.

$Y_n$ $1/Y_n$ $n \rightarrow \infty$ $Y_n$ $Y_n = 0$ $e$

— Cliff AB
źródło

2

e

$e$

1

@ whuber: Nie korzystałem z Box-Mullera z powodu wymaganej transformacji dziennika zbyt bezpośrednio, aby wykładniczy w mojej książce. I to już odruchowo dozwolone cos i grzechu, ale to było tylko dlatego, że zapomniał o złożonej analizy przez chwilę, więc dobry punkt.

— Cliff AB,

1

n_{1}

$n_1$

n_{2}

$n_2$

ϕ (\sqrt{2})

$\phi(\sqrt{2})$

n_{1}

$n_1$

n_{2}

$n_2$

\infty

$\infty$

e

$e$

n_{2}

$n_2$

n_{1}

$n_1$

2

Oto inny sposób, w jaki można to zrobić, choć jest on dość wolny. Nie rości sobie pretensji do skuteczności, ale oferuję tę alternatywę w duchu kompletności.

$n$ $U_1, \cdots , U_n \sim \text{IID U}(0,1)$ $e$ $^\dagger$

mi (\frac{ja (U_{ja} ⩾ 1 / mi)}{U_{ja}}) = \int_{1 / mi}^{1} \frac{re u}{u} = 1.

$\mathbb{E} \Bigg( \frac{\mathbb{I}(U_i \geqslant 1 / e) }{U_i} \Bigg) = \int \limits_{1/e}^1 \frac{du}{u} = 1.$

$e$ $u_{(1)} \geqslant \cdots \geqslant u_{(n)}$

S_{n} (k) \equiv \frac{1}{n} \sum_{i = 1}^{k} \frac{1}{u_{(i)}} for all k = 1, . ., n .

$S_n(k) \equiv \frac{1}{n} \sum_{i=1}^k \frac{1}{u_{(i)}} \quad \text{for all } k = 1, .., n.$

$m \equiv \min \{ k | S(k) \geqslant 1 \}$ $1/e$ $e$

\hat{e} \equiv \frac{2}{u_{(m)} + u_{(m + 1)}} .

$\hat{e} \equiv \frac{2}{u_{(m)} + u_{(m+1)}}.$

$1/e$ $e$ $e$

Implementacja w R: Metodę można zaimplementować przy Rużyciu runifdo generowania jednolitych wartości. Kod jest następujący:

EST_EULER <- function(n) { U <- sort(runif(n), decreasing = TRUE);
                           S <- cumsum(1/U)/n;
                           m <- min(which(S >= 1));
                           2/(U[m-1]+U[m]); }

$e$

set.seed(1234);

EST_EULER(10^3);
[1] 2.715426

EST_EULER(10^4);
[1] 2.678373

EST_EULER(10^5);
[1] 2.722868

EST_EULER(10^6); 
[1] 2.722207

EST_EULER(10^7);
[1] 2.718775

EST_EULER(10^8);
[1] 2.718434

> exp(1)
[1] 2.718282

$^\dagger$ $e$

— Przywróć Monikę
źródło

Przybliżone za pomocą symulacji Monte Carlo

Oszacowanie przez ładowanie początkoweeee

Wyniki

Oszacowanie przez ładowanie początkowe $e$