Il s'agit d'un problème célèbre connu sous le nom de Gull's Lighthouse, d'après un exemple de Gull en 1988. Il a de profondes implications lorsqu'on franchit une étape supplémentaire dans les sciences sociales et la physique. Vous avez en fait suffisamment d'informations pour résoudre ce problème grâce à des tests d'acceptation-rejet, mais si vous souhaitez utiliser MCMC, n'hésitez pas.

Examinons votre problème d'une autre manière, sachant que le phare est à 100 m du rivage est en fait beaucoup d'informations. Ce morceau est un gros problème, sinon vous auriez à résoudre à la fois la distance et l'emplacement.

En remarque, la distribution de Cauchy est large de 636 plages semi-interquartiles à 99,99% HDR, donc si vous cherchiez intervalles étroits, vous pouvez presque l'oublier. Cependant, si vous utilisez une approximation normale, elle ne convergera jamais car ce problème viole le théorème de la limite centrale et si vous estimez $ \ sigma $ comme écart type, alors si $ n $ est la taille de l'échantillon, alors $ s ^ 2 \ à \ infty $ comme $ n \ à \ infty $. Si vous faisiez cela, vos estimations deviendraient de plus en plus mauvaises à mesure que vous ajoutiez des dates.

Vous devez noter que vous observez l'ensemble $ \ {x_i, i \ in {1 \ dots {n} } \} $ et il y a un point $ \ mu $ où le phare serait perpendiculaire à la plage. Notez que ce n'est pas le positionnement des $ x_i $ qui sont uniformément distribués, mais plutôt l'angle, $ \ theta $. Nous n'avons pas observé l'angle, seulement l'ensemble des points.

La relation des points à l'angle est $$ \ frac {x_i- \ mu} {100} = \ tan \ theta. $$ En termes d'angle, cela peut être exprimé par $$ \ theta = \ tan ^ {- 1} \ left (\ frac {x_i- \ mu} {100} \ right) $$

Depuis $$ p (x_i | \ mu, 100) = p (\ theta | \ mu, 100) | \ frac {\ mathrm {d} \ theta} {\ mathrm {d} x_i} | $$ et $ \ theta $ est uniforme sur la région $ (- \ pi / 2, \ pi / 2) $ on sait que $$ p (x_i | \ mu, 100) = \ frac {1} {\ pi} | \ frac {\ mathrm { d} \ theta} {\ mathrm {d} x_i} |, $$ qui vaut $$ p (x_i | \ mu, 100) = \ frac {1} {\ pi} \ frac {100} {100 ^ 2 + (x_i- \ mu) ^ 2}. $$

Donc, d'après le théorème de Bayes, si le point de rivage perpendiculaire au phare est $ (\ mu, 0) $, alors le phare est situé à $ (\ mu, 100) $. L'estimateur de ce point, puisque vous n'avez aucune information préalable sur $ \ mu $ est $$ p (\ mu | x_1 \ dots {x_n}, 100) \ propto \ prod_ {i = 1} ^ n \ frac {1} {\ pi} \ frac {100} {100 ^ 2 + (x_i- \ mu) ^ 2}, \ forall \ mu \ in \ Re. $$

Si vous voulez rendre le problème plus intéressant, supposons que vous ne savez pas que l'emplacement est à cent mètres du rivage. Votre problème devient alors le problème bidimensionnel le plus adapté à stan. Votre problème devient alors $$ p (\ mu, \ gamma | x_1 \ dots {x_n}) \ propto \ prod_ {i = 1} ^ n \ frac {1} {\ pi} \ frac {\ gamma} {\ gamma ^ 2 + (x_i- \ mu) ^ 2}, \ forall \ mu \ in \ Re; \ gamma \ in \ Re ^ {++}. $$

En guise de remarque, ce problème est liée à la célèbre sorcière d'Agnèsi, étudiée par Maria Agnesi et publiée en 1748 et au problème plus fondamental étudié par Fermat et Grandi en 1703. Elle l'avait inclus dans le tout premier véritable manuel de calcul. Il a fourni aux étudiants une éducation de l'algèbre à travers des équations intégrales et différentielles.

Les positions des flashs sont modélisées avec précision par le Cauchy, comme l'a dit @ user25459, donc ce n'est pas un hasard si Stan est capable d'estimer les vraies valeurs en utilisant ce modèle (j'utilise $ \ alpha $ pour $ \ mu $, et $ \ beta $ pour $ \ gamma $):

  model <- 'data {int<lower = 0> N; réel x_ [N]; } paramètres {alpha réel; real<lower = 0> beta; } modèle {alpha ~ uniforme (0, 20); beta ~ uniforme (0, 50); pour (k en 1: N) {x_ [k] ~ cauchy (alpha, bêta); }} ' 

Essayons ceci dans Stan:

  alpha <- 10 # unknown true valuesbeta <- 30 ######### ######### set.seed (123) N <- 100theta_k <- runif (N, -pi / 2, pi / 2) x_k <- beta * tan (theta_k) + alphastan_input <- liste (x_ = x_k, N = N) fit <- stan (model_code = model, data = stan_input, iter = 1000, verbose = FALSE) fit2 <- stan (fit = fit, data = stan_input, iter = 5000, warmup = 2000, verbose = FALSE) print (fit2, pars = c ("alpha", "beta"))  

Obtenez-nous:

  Inférence pour le modèle Stan: 36eef3c3637fb3a9564529926f8463fe .4 chaînes, chacune avec iter = 5000; échauffement = 2000; mince = 1; tirages post-échauffement par chaîne = 3000, nombre total de tirages post-échauffement = 12000. moyenne se_mean sd 2,5% 25% 50% 75% 97,5% n_eff Rhatalpha 8,96 0,05 3,91 1,61 6,13 8,95 11,68 16,80 7079 1bêta 30,29 0,05 4,18 22,84 27,45 29,98 32,88 39,18 8470 1 Les échantillons ont été prélevés à l'aide de NUTS (diag_e) le mercredi 11 janvier 18:28: 57 2017.Pour chaque paramètre, n_eff est une mesure brute de la taille effective de l'échantillon et Rhat est le facteur de réduction d'échelle potentiel sur les chaînes fractionnées (à la convergence, Rhat = 1).  

l'un des exemples:

  la <- extract (fit2) alpha_hats <- as.numeric (la $ alpha) # la ligne verticale est la valeur vraie  

À propos des angles, notez que

$$ \ theta_k = \ arctan \ frac {x_k - \ alpha} {\ beta} $ $

Chaque $ \ theta_k $ est une fonction de deux variables aléatoires, donc est elle-même une variable aléatoire.

Pour trouver la distribution empirique de l'angle $ \ theta_k $, nous utilisons les échantillons aléatoires de $ \ alpha$ et $ \ beta $.

  theta_1 <atan ((x_k [1] - alpha_hats) / beta_hats) hist (theta_1, breaks = 100, xlab = bquote (theta [1]), ylab = "", main = "Distribution pour le 1er angle", yaxt = 'n')