密度関数が $1/(1+x^2)$ に比例する乱数（独立である必要はない）を10万個作る問題です。『Rで楽しむベイズ統計入門』に載っているRのプログラムは（ここで使わない採択を数えるカウンタは省略して）次のようになっていました：
N = 100000 # 繰返し回数
a = numeric(N) # 値を保存するための長さNの配列
x = 0 # 初期値
p = 1 / (1 + x^2) # 確率
for (i in 1:N) {
y = rnorm(1, x, 1) # 候補（選び方は対称なら何でもよい）
q = 1 / (1 + y^2) # 確率
if (runif(1) < q/p) { # 更新
x = y
p = q
}
a[i] = x
}
実際にこれをRにコピペしてやってみるとわかるように，今どきのコンピュータなら10万個の計算はほぼ一瞬で終わります。
得られた乱数のヒストグラムを描いてみましょう：
hist(a, freq=FALSE, breaks=seq(-1000,1000,0.2), xlim=c(-5,5), col="gray")
コーシー分布の密度関数を重ね書きしてみましょう：
curve(dcauchy(x), add=TRUE)
ぴったり重なるはずです。結果は，本書の図7.3右のようになります。
同じことをStanでやってみましょう。
まず，次のようなStanコードを作り，例えば ex72.stan という名前で，Rの作業ディレクトリに保存します。
parameters {
real x;
}
model {
target += -log(1 + square(x));
}
テキストエディタで最後の } を打ち込んだ後に Enter キーを打って，カーソルが } の次の行の先頭に行く状態にして，保存してください。
改行で終わらない行は incomplete line と呼ばれ，プログラミングの世界では嫌われます。Stan でもときどき make sure Stan code ends with a blank line というメッセージが出ます。
上のコードは，パラメータとして実数（real は C++ の double に相当）の x を使うことと，モデルとして確率密度の対数（定数項を除く）が $-\log(1+x^2)$ であることを表します。target += は確率密度の対数を追加する命令です。決まり文句ですので，target = target - log(1 + square(x)) とか target -= log(1 + square(x)) とかに書き直してはいけません。square(x) は x * x や x^2 と同じことです。
そして，例えば次のようにRに打ち込みます：
fit = stan("ex72.stan")
この stan() という関数は，デフォルトでは4本のシミュレーションを（可能なら並行して）実行します（chains=4）。1本あたり，デフォルトでは2000回繰り返しますが，その半分はウォームアップ（バーンイン）です。つまり，1本あたり1000個のサンプル値を作りますので，全部で4000個が得られます。
これを打ち込むと，C++への変換とコンパイルでかなり長い無反応な時間が続きますが，その後で計算の経過が出力されます。最後に警告メッセージも出るかもしれませんが，とりあえず無視します。
rstan_options(auto_write=TRUE) にしている場合，ex72.stan をコンパイルした結果は ex72.rds というファイルに入ります。ex72.stan（のハッシュ値）が変わったら，コンパイルし直してくれるはずです。
fit に代入された結果を表示します：
fit
Inference for Stan model: ex72.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
x -0.66 0.47 12.54 -11.53 -1.11 -0.05 0.81 8.57 722 1.00
lp__ -1.25 0.08 1.60 -5.96 -1.64 -0.81 -0.17 0.00 379 1.01
Samples were drawn using NUTS(diag_e) at Sun Feb 4 17:32:10 2018.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
x（と lp__ つまり対数事後分布 log posterior）について，平均とその標準誤差，標準偏差，％点，有効サンプルサイズ $n_{\text{eff}}$，それに $\hat{R}$ という統計量が出力されます。これらの詳しい定義は Bayesian Data Analysis, 3rd edition の11.4〜11.5に書かれています。$\hat{R} \approx 1$ なら収束したと考えられます。$n_{\text{eff}}$ が小さければ，そのパラメータの推定誤差が大きいので，後述のようにサンプルサイズを増します。
警告メッセージについては，説明へのリンクがありますので，そちらを見て，対処できるものは対処すればよいでしょう。例えば There were XXX transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded と出たら，リンクをたどって問題点を見つけ，
fit = stan("ex72.stan", control=list(max_treedepth=20))
のようにオプションでパラメータを変更します。ただ，うまく収束しなかったり警告が出たりすることは，必ずしもパラメータの変更で解消できるわけではありません。
とりあえずうまくいくようなら，サンプルサイズを大きくしてみましょう。このページの最初に挙げたRコードでは10万でしたので，それに合わせるために，シミュレーション1本あたりの（ウォームアップを含む）繰返し数 iter とウォームアップ（バーンイン）の数 warmup を調整します。例えば
fit = stan("ex72.stan", iter=26000, warmup=1000, control=list(max_treedepth=20))
とすれば，1本あたりのサンプルサイズが 25000 になりますので，4本でちょうど10万個が得られます。
より詳しく調べるために，4本のシミュレーション（ウォームアップを除く）の結果をつなぎ合わせた100000個の x や lp__ の値を取り出しましょう：
e = extract(fit)
これで e$x に100000個の x の値が取り出されます。
本書と同様にヒストグラム表示してみましょう。とてつもなく絶対値の大きな値が混じっていることがあるので，どんな場合にもうまくいくように，少し工夫しています。
x = ifelse(abs(e$x) < 100, e$x, 100) # 絶対値を100以内に止める hist(x, freq=FALSE, breaks=seq(-100,100,0.2), xlim=c(-5,5), col="gray") curve(dcauchy(x), add=TRUE)
結果のグラフは本書のRコードによる結果（図7.3右）と同じです。コーシー分布ですから，ときどきとてつもなく絶対値の大きな値が混じります。この部分が，本書のRコードより上のStanの結果のほうが正確です。これは，本書のRコードで次の候補を選ぶ y = rnorm(1, x, 1) の跳び幅が狭く，あまり絶対値の大きな値が出にくいためです。
extract() で取り出したサンプルは，デフォルトではランダムにシャフルされています。シャフルしないそのままの順番で4本のシミュレーション結果をつないだものを得るには，extract() に permuted=FALSE というオプションを付けるか，あるいは
e = as.data.frame(fit)
とします。ウォームアップを含むまったくの生の値は fit@sim$samples[[1]]$x から fit@sim$samples[[4]]$x までに入っています（chains=4 の場合）。
関数 stan() の実行は，最初はコンパイルの時間が長くかかりますが，次からはずっと速くなります。ただ，stan() を別の関数の中で使うと，毎回コンパイルが起こってしまい，遅くなります。その場合は例えば
model = stan_model("ex72.stan")
f = function(data) {
fit = sampling(model, iter=26000, warmup=1000)
mean(extract(fit)$x)
}
のように stan_model() と sampling() に分けて計算します。