39  Intervalos de Confiança - Aprofundamento

Seja \(\boldsymbol{X}_n = (X_1,\dots,X_n)\) amostra aleatória e \(X \sim f_\theta, \theta \in \Theta\).

Dizemos que \([I_1(\boldsymbol{X}_n), I_2(\boldsymbol{X}_n)]\) é um intervalo de confiança exato para \(g(\theta)\) com coeficiente de confiança \(\gamma \in (0,1)\) se, e somente se: \[ P_\theta\left(I_1(\boldsymbol{X}_n) \leq g(\theta) \leq I_2(\boldsymbol{X}_n)\right) = \gamma, \forall \theta \in \Theta. \] em que \(I_1(\boldsymbol{X}_n), I_2(\boldsymbol{X}_n)\) são estatísticas.

Observação

Se \(P_\theta\left(I_1(\boldsymbol{X}_n) \leq g(\theta) \leq I_2(\boldsymbol{X}_n)\right) \geq \gamma, \forall \theta \in \Theta\), então \([I_1(\boldsymbol{X}_n), I_2(\boldsymbol{X}_n)]\) é um Intervalo de Confiança (IC) de pelo menos \(\gamma\)

Notação

\[ \mathrm{IC}(g(\theta), \gamma) = [I_1(\boldsymbol{X}_n), I_2(\boldsymbol{X}_n)] \]

Observe que \(\mathrm{IC}(g(\theta),\gamma)\) é um intervalo aleatório que não depende de “\(\theta\)”.

Na prática, observamos a amostra \(\boldsymbol{x}_n = (x_1, \dots, x_n)\) e calculamos o IC observado

Notação

\[ \mathrm{IC}_{\mathrm{Obs}}(g(\theta), \gamma) = [I_1(\boldsymbol{x}_n), I_2(\boldsymbol{x}_n)] \]

Observe que este é um intervalo numérico

Portanto, \[ P_\theta\left(I_1(\boldsymbol{X}_n) \leq g(\theta) \leq I_2(\boldsymbol{X}_n)\right) = \begin{cases} 1, g(\theta) \in \mathrm{IC}_{\mathrm{Obs}}(g(\theta), \gamma) \\ 0, \mathrm{c.c.} \end{cases} \]

39.1 Interpretação em termos de repetições:

Se repetirmos o experimento, mantendo as mesmas condições, então esperamos que em \(\gamma \cdot 100\%\) dos experimentos os ICs contenham \(g(\theta)\). Em outras palavras, \[ \begin{aligned} \#\frac{(g(\theta) \in \mathrm{IC}_{\mathrm{Obs}})}{N} \approx \gamma \\ \left[ \frac{1}{N} \sum \mathbb{1}_{\{\mathrm{IC}_{\mathrm{Obs}}^{(i)}\}}(g(\theta)) \stackrel{N \uparrow \infty}{\rightarrow} \gamma \right] \end{aligned} \]

Dizemos que \[ \frac{1}{N} \sum \mathbb{1}_{\{\mathrm{IC}_{\mathrm{Obs}}^{(i)}\}}(g(\theta)) \] é a cobertura de \(\mathrm{IC}(g(\theta),\gamma)\)

O que dizer sobre \(\mathrm{IC}_{\mathrm{Obs}}\) em relação a \(g(\theta)\)?

Temos uma confiança de \(\gamma\cdot 100\%\) de que \(g(\theta) \in \mathrm{IC}_{\mathrm{Obs}}\).

39.2 Exemplos

39.2.1 Exemplo Normal

Seja \(\boldsymbol{X}_n\) a.a. de \(X \sim N(\mu,\sigma^2), \theta = (\mu, \sigma^2) \in \Theta = \mathbb{R}\times\mathbb{R}^+\). Encontre um IC com coeficiente de confiança \(\gamma = 95\%\)…

…para \(g(\theta = \mu)\)

Sabemos que \[ \bar{X} = \frac{1}{n} \sum X_i \sim N(\mu, \sigma^2/n). \] Além disso, \[ \sum \frac{(X_i - \bar{X})^2}{\sigma^2} \sim \chi^2_{n-1} \]

Por definição de t-student, com as V.As das distribuições independentes, \[ \begin{aligned} \frac{N(0,1)}{\sqrt{\frac{\chi^2_k}{k}}} &\sim t_k \\ \Rightarrow\frac{\sqrt{n} \frac{(\bar{X} - \mu)}{\sigma}}{\sqrt{\frac{\sum \frac{(X_i - \bar{X})^2}{\sigma^2}}{n-1}}} &\sim t_{(n-1)} \\ \Rightarrow \frac{\sqrt{n} (\bar{X} - \mu)}{\sqrt{S^2_{n-1}(\boldsymbol{X}_n)}} &\sim t_{(n-1)} \end{aligned} \] logo, podemos sempre encontrar \(c_{1, \gamma}, c_{2, \gamma}\) tais que \[ P_\theta\left(c_{2, \gamma} \leq \frac{\sqrt{n} (\bar{X} - \mu)}{\sqrt{S^2_{n-1}(\boldsymbol{X}_n)}} \leq c_{1, \gamma}\right) = \gamma \]

Note que \[ c_{2, \gamma} \leq \frac{\sqrt{n} (\bar{X} - \mu)}{\sqrt{S^2_{n-1}(\boldsymbol{X}_n)}} \leq c_{1, \gamma} \iff \bar{X} - c_{2,\gamma} \sqrt{\frac{S^2_{n-1}(\boldsymbol{X}_n)}{n}} \leq \mu \leq \bar{X} - c_{1,\gamma} \sqrt{\frac{S^2_{n-1}(\boldsymbol{X}_n)}{n}} \]

Portanto, \[ P_\theta\left(\bar{X} - c_{2,\gamma} \sqrt{\frac{S^2_{n-1}(\boldsymbol{X}_n)}{n}} \leq g(\theta) \leq \bar{X} - c_{1,\gamma} \sqrt{\frac{S^2_{n-1}(\boldsymbol{X}_n)}{n}}\right) = \gamma \]

Pela definição de Intervalo de Confiança, \[ \mathrm{IC}(\mu,\gamma) = \left[ \bar{X} - c_{2,\gamma} \sqrt{\frac{S^2_{n-1}(\boldsymbol{X}_n)}{n}}, \bar{X} - c_{1,\gamma} \sqrt{\frac{S^2_{n-1}(\boldsymbol{X}_n)}{n}} \right] \] em que \[ S^2_{n-1}(\boldsymbol{X}_n) = \frac{1}{n-1} \sum (X_i - \bar{X})^2 \] e \(c_{1,\gamma}, c_{2,\gamma}\) são os quantis obtidos da distribuição t-student com \(n-1\) graus de liberdade que satisfaçam \[ P_\theta\left(\bar{X} - c_{2,\gamma} \sqrt{\frac{S^2_{n-1}(\boldsymbol{X}_n)}{n}} \leq g(\theta) \leq \bar{X} - c_{1,\gamma} \sqrt{\frac{S^2_{n-1}(\boldsymbol{X}_n)}{n}}\right) = \gamma \] no caso simétrico (minimiza o IC para distribuições como a Normal), \(c_{2,\gamma} = - c_{1,\gamma}\). Note que não é possível construir ICs simétricos dessa forma para distribuições estritamente positivas, como a qui-quadrado.

39.3 Quantidades Pivotais

Dizemos que \(Q(g(\theta), \boldsymbol{X}_n)\) é uma quantidade pivotal para \(g(\theta)\) se, e somente se,

  1. \(Q(g(\theta), \boldsymbol{X}_n)\) depende de \(g(\theta)\)

  2. A distribuição de \(Q(g(\theta), \boldsymbol{X}_n)\) não depende de “\(\theta\)”

  3. Existem \(a_1, a_2\), que não dependem de \(g(\theta)\), tais que \(c_1 \leq Q(g(\theta),\boldsymbol{X}_n) \leq c_2 \iff a_1 \leq g(\theta) \leq a_2\)

39.3.1 Exemplos

\(X \sim N(\mu, \sigma^2)\)

  1. Se \(\sigma^2\) é conhecido e \(\theta = \mu. g(\theta) = \mu\), então uma quantidade pivotal é dada por \[ Q(\mu, \boldsymbol{X}_n) = \sqrt{n}\frac{\bar{X} - \mu}{\sqrt{\sigma^2}} \sim N(0,1) \]

\[ \mathrm{IC}(\mu,\gamma) = \bar{X} \mp c_{\gamma} \sqrt{\frac{\sigma^2}{n}} \]

  1. Se \(\sigma^2\) é desconhecido e \(\theta = (\mu, \sigma^2), g(\theta) = \mu\), então uma quantidade pivotal é dada por \[ Q(\mu, \boldsymbol{X}_n) = \sqrt{n}\frac{\bar{X}-\mu}{\sqrt{S^2_{n-1}(\boldsymbol{X}_n)}} \sim t_{(n-1)} \]

\[ \mathrm{IC}(\mu,\gamma) = \bar{X} \mp c_{\gamma} \sqrt{\frac{S^2_{n-1}(\boldsymbol{X}_n)}{n}} \] essa é também uma quantidade pivotal para 1., mas o contrário não vale.

  1. Se \(\mu\) é conhecido e \(\theta = \sigma^2, g(\theta) = \sigma^2\), então

\[ Q(\sigma^2, \boldsymbol{X}_n) = \sum \frac{(X_i - \bar{X})^2}{\sigma^2} \sim \chi^2_{n} \]

\[ \mathrm{IC}(\sigma^2,\gamma) = \left[\frac{\sum (X_i - \mu)^2}{c_{2,\gamma}}, \frac{\sum (X_i - \mu)^2}{c_{1,\gamma}}\right] \]

  1. Se \(\mu\) é desconhecido e \(\theta = (\mu, \sigma^2), g(\theta) = \sigma^2\), então \[ Q(\sigma^2, \boldsymbol{X}_n) = \sum \frac{(X_i - \bar{X})}{\sigma^2} \sim \chi^2_{(n-1)} \]

\[ \mathrm{IC}(\sigma^2,\gamma) = \left[\frac{\sum (X_i - \bar{X})^2}{c_{2,\gamma}}, \frac{\sum (X_i - \bar{X})^2}{c_{1,\gamma}}\right] \]

39.3.1.1 Exponencial

Seja \(X \sim \mathrm{Exp}(\theta), \theta > 0\). Encontre uma quantidade pivotal para \(\theta\).

Note que \(\sum X_i \sim \mathrm{Gama}(n,\theta)\) com F.G.M. dada por \[ M_{\sum X_i}(t) = \left(\frac{\theta}{\theta-t}\right)^n \] note ainda que \[ \begin{aligned} M_{\sum X_i}(t) &= E(\mathrm{e}^{t \sum X_i}) = \left(\frac{\theta}{\theta-t}\right)^n \\ \Rightarrow M_{\theta \sum X_i}(t) &= E(\mathrm{e}^{t\theta \sum X_i}) \\ &=M_{\sum X_i}(t \theta) = \left(\frac{\theta}{\theta-\theta t}\right)^n = \left(\frac{1}{1-t}\right)^n \\ \Rightarrow M_{2\theta\sum X_i}(t) &= \left(\frac{1}{1-2t}\right)^n \\ \Rightarrow &2\theta \sum X_i \sim \chi^2_{(2n)} \end{aligned} \]

Portanto, \(Q(\theta, \boldsymbol{X}_n) = 2\theta \sum X_i\) é uma quantidade de interesse para \(\theta\)

\[ \begin{aligned} c_{1,\gamma} &\leq Q(\theta, \boldsymbol{X}_n) \leq c_{2,\gamma} \\ \iff c_{1,\gamma} &\leq 2\theta \sum X_i \leq c_{2\gamma} \iff \frac{c_{1,\gamma}}{2\sum X_i} \leq \theta \leq \frac{c_{2,\gamma}}{2\sum X_i} \end{aligned} \]

39.3.1.2 Uniforme

Seja \(X \sim \mathrm{Unif}(0, \theta), \theta > 0\). Encontre uma quantidade pivotal para \(\theta\).

Note que \(X_{(n)} = \max \boldsymbol{X}_n\) é uma estatística suficiente cuja f.d.p. é dada por \[ f_(\theta)^{\boldsymbol{X}_{(n)}}(x) = \frac{n x^{n-1}}{\theta^n} \mathrm{1}_{(0,\theta]}(x) \]

Seja \(Y = \frac{X_{(n)}}{\theta}\), então \[ f_{\theta}^{Y}(y) = f_{\theta}^{\boldsymbol{X}_{(n)}}(y \theta) \cdot |J| \] em que \(J = \theta\) (determinante jacobiano)

\[ f_{\theta}^{Y}(y) = \frac{n (y\theta)^{n-1}}{\theta^n} \theta \mathbb{1}_{(0,\theta]}(y \theta) = n y^{n-1} \mathbb{1}_{(0,1]}(y) \] não depende de “\(\theta\)”! Logo, \(Q(\theta, \boldsymbol{X}_n) = \frac{\boldsymbol{X}_{(n)}}{\theta}\) é uma quantidade pivotal para \(\theta\).

\[ \begin{aligned} \frac{\boldsymbol{X}_{(n)}}{c_{2,\gamma}} \leq \theta \leq \frac{\boldsymbol{X}_{(n)}}{c_{1,\gamma}} \Rightarrow \mathrm{IC}(\theta,\gamma) = \left[ \frac{\boldsymbol{X}_{(n)}}{c_{2,\gamma}}, \frac{\boldsymbol{X}_{(n)}}{c_{1,\gamma}} \right] \end{aligned} \] em que os quantis são obtidos da distribuição de \(Y\), nesse caso, \(Y \sim \mathrm{Beta}(n,1)\):

\[ \begin{aligned} \int^{c_{1,\gamma}}_0 ny^{n-1} dy = \frac{1-\gamma}{2}&\ \ \ \int_{c_{2,\gamma}}^{1} ny^{n-1}dy = \frac{1-\gamma}{2} \\ \Rightarrow \begin{cases} y^n \rvert_0^{c_{1,\gamma}} = \frac{1-\gamma}{2} \\ y^n \rvert^1_{c_{2,\gamma}} = \frac{1-\gamma}{2} \\ \end{cases} \Rightarrow c_{1,\gamma} &= \left(\frac{1-\gamma}{2}\right)^{1/2}\ \ \ \ c_{2,\gamma} = \left(\frac{1+\gamma}{2}\right)^{1/2} \end{aligned} \]

Quantis que minimizam a amplitude do IC

Seja \(Q(g(\theta), \boldsymbol{X}_n)\) uma quantidade pivotal com função densidade de probabilidade \(f\).

\(c_{1,\gamma}, c_{2,\gamma}\) devem ser obtidos \[ \int_{c_{1,\gamma}}^{c_{2,\gamma}} f(x) dx = \gamma \tag{39.1}\]

Note que em geral infinitas combinações desses quantis satisfazem (39.1). Podemos usar o par que satisfaz \[ \int^{c_{1,\gamma}}_{-\infty} f(y) dy = \frac{1-\gamma}{2}\ \ \mathrm{e} \ \ \int_{c_{2,\gamma}}^{\infty} f(y)dy = \frac{1-\gamma}{2} \]

Esse método produz um intervalo de confiança simétrico, não necessariamente o de menor amplitude, mas é mais fácil de encontrar. O intervalo de confiança com menor amplitude com quantis que satisfaçam (39.1) é obtido minimizando \(|c_{2,\gamma} - c_{1,\gamma}|\) sujeito a (39.1).

Se \(f\) for unimodal e bicaudal, pode-se demonstrar que \(c_{1,\gamma}, c_2{\gamma}\) que produzem amplitude mínima e satisfazem (39.1) são tais que \[ f(c_{1,\gamma}) = f(c_{2,\gamma}) \] ou seja, tem mesma densidade.

using Distributions, Random, StatsBase, LaTeXStrings

Random.seed!(8)

theta0 = 4
n = 10
MC = 10000

I1 = []
I2 = []
for _ in 1:MC
    d = Exponential(1 / theta0) # Parâmetro média => 1/theta parâmetro taxa
    x = rand(d, n)

    # Construindo um IC com 95% de confiança
    gamma = 0.95
    quiquadrado = Chisq(2 * n)
    c1 = quantile(quiquadrado, (1 - gamma) / 2)
    c2 = quantile(quiquadrado, gamma + (1 - gamma) / 2)
    push!(I1, round(c1 / (2 * sum(x)), digits = 4))
    push!(I2, round(c2 / (2 * sum(x)), digits = 4))
end

acertos = [I1 .<= theta0 .<= I2]
display(L"\mathrm{IC}_{\mathrm{Obs (1)}}(\theta,0.95)=[%$(I1[1]), %$(I2[1])]")
display(L"\text{ICs que contém}\ \theta: %$((sum(acertos[1]))/MC * 100)\%")

\(\mathrm{IC}_{\mathrm{Obs (1)}}(\theta,0.95)=[2.7345, 9.7424]\)

\(\text{ICs que contém}\ \theta: 94.97\%\)

39.3.2 Quantidades pivotais aproximadas ou assintóticas

Dizemos que \(Q(g(\theta), \boldsymbol{X}_n)\) é uma quantidade pivotal aproximada ou assintótica se, e somente se

  1. \(Q(g(\theta), \boldsymbol{X}_n)\) depende de \(g(\theta)\);

  2. A distribuição assintótica de \(Q(g(\theta), \boldsymbol{X}_n)\) não depende de \(g(\theta)\);

  3. Existem \(a_1, a_2\), que não dependem de \(g(\theta)\), tais que \(c_1 \leq Q(g(\theta),\boldsymbol{X}_n) \leq c_2 \iff a_1 \leq g(\theta) \leq a_2\). Esses valores dependem apenas de \(\boldsymbol{X}_n\).

Podemos usar o TLC para estimadores:

Se \(T(\boldsymbol{X}_n)\) for um estimador para \(\theta\) assintoticamente normal, então

\[ \sqrt{n} (T(\boldsymbol{X}_n) - \theta) \stackrel{\mathcal{D}}{\rightarrow} N_p(0, V_\theta), \forall \theta \in \Theta. \] em que \(V_\theta\) é uma matriz positiva definida. Para \(p=1\), \[ \sqrt{n} (T(\boldsymbol{X}_n) - \theta) \stackrel{\mathcal{D}}{\rightarrow} N_1(0, V_\theta), \forall \theta \in \Theta. \] em que \(V_\theta > 0\).

Se \(g: \Theta \rightarrow \mathbb{R}\) for uma função tal que \(g'(\theta) \neq 0\) \(g'(\theta)\) é contínua, então, \[ \sqrt{n} (g(T(\boldsymbol{X}_n)) - g(\theta)) \stackrel{\mathcal{D}}{\rightarrow} N_1(0, g'(\theta)^2V_\theta), \forall \theta \in \Theta. \]

Além disso,

  1. \[ \frac{\sqrt{n} (g(T(\boldsymbol{X}_n)) - g(\theta))}{\sqrt{g'(\theta)^{2}V_\theta}} \stackrel{\mathcal{D}}{\rightarrow} N_1(0, 1), \forall \theta \in \Theta. \]

  2. Pelo teorema de Slutsky \[ \frac{\sqrt{n} (g(T(\boldsymbol{X}_n)) - g(\theta))}{\sqrt{g'(T(\boldsymbol{X}_n))^{2}V_{T(\boldsymbol{X}_n)}}} \stackrel{\mathcal{D}}{\rightarrow} N_1(0, 1), \forall \theta \in \Theta. \]

Ou seja, \[ Q(g(\theta), \boldsymbol{X}_n) = \frac{\sqrt{n} (g(T(\boldsymbol{X}_n)) - g(\theta))}{\sqrt{g'(T(\boldsymbol{X}_n))^{2}V_{T(\boldsymbol{X}_n)}}} \] é uma quantidade pivotal aproximada para \(g(\theta)\). Note que, com \(T(\boldsymbol{X}_n) = \hat(\theta)\), \[ \begin{aligned} &c_1 \leq Q(g(\theta), \boldsymbol{X}_n) \leq c_2 \\ \iff& c_1 \leq \frac{\sqrt{n} (g(T(\boldsymbol{X}_n)) - g(\theta))}{\sqrt{g'(T(\boldsymbol{X}_n))^{2}V_{T(\boldsymbol{X}_n)}}} \leq c_2 \\ \iff& c_1 \sqrt{\frac{g'(\hat\theta)^2}{n} V_{\hat\theta}} \leq g(\hat\theta) - g(\theta) \leq c_2\sqrt{\frac{g'(\hat\theta)^2}{n} V_{\hat\theta}} \\ \iff& g(\hat\theta) - c_2 \sqrt{\frac{g'(\hat\theta)^2}{n} V_{\hat\theta}} \leq g(\theta) \leq g(\hat\theta) - c_1\sqrt{\frac{g'(\hat\theta)^2}{n} V_{\hat\theta}} \\ \end{aligned} \] Tomando \(c_1 = -c_2\), pois a distribuição assintótica é normal, temos \[ \mathrm{IC}^a(g(\theta),\gamma) = \left[ g(\hat\theta) - c_2 \sqrt{\frac{g'(\hat\theta)^2}{n} V_{\hat\theta}}, g(\hat\theta) + c_2 \sqrt{\frac{g'(\hat\theta)^2}{n} V_{\hat\theta}} \right] \] em que \(c_2\) é obtido dos quantis da normal padrão tal que \[ P(-c_2 \leq N(0,1) \leq c_2) = \gamma \]

Erro padrão

A quantidade \[ \sqrt{\frac{g'(\hat\theta)^2}{n} V_{\hat\theta}} \] é o erro-padrão do estimador (seu desvio padrão).

Vamos conferir por simulações de monte carlo e o método de Newton-Raphson

  1. Para \(\theta\)
  2. Para \(g(\theta) = P(X < 950) \Rightarrow g'(\theta) = f_\theta(950)\)
  3. Para \(g(\theta) = \ln f_\theta(x) g'(\theta) = - \psi_1(\theta) + \ln(\theta) - 1/\theta\)
  4. Para \(g(\theta) = f_\theta(x)\)
using SpecialFunctions, Distributions, Random, LaTeXStrings
# Modelo Gama(θ, 1)
# IC assintótico
# Calcule o IC
#a-) para θ
#b-) para g(θ) = P(X < 950) => g'(θ) = f_θ(950)
#c-) para g_x(θ) = ln f_θ(x) g'(θ) = - ψ_1(θ) + ln(θ) - 1/θ. x = θ
#d-) para g_x(θ) = f_\theta(x). x = θ


# Sabemos que θ.MV ~a~ N(θ, I_n(θ)^-1)
# => IC^a(θ,γ) = θ.MV ∓ c_2 * √(I_1(θ))

Random.seed!(13)

function newton_raphson(x)

    theta::Vector{Float64} = []
    append!(theta, mean(x)) # Chute inicial = média
    erromax = 10^(-5)
    erro = Inf
    i = 1
    # iteracoesMax = 6 # Podemos também definir apenas um erro máximo
    while erro > erromax # && i < iteracoesMax
        append!(
            theta, theta[i] - (sum(log.(x)) - n * digamma(theta[i])) /
                (-n * trigamma(theta[i]))
        )
        erro = abs(theta[i + 1] - theta[i])
        # println("Erro na iteração $i: $erro")
        i += 1
    end
    # println("Theta final: $(theta[length(theta)])")
    # println("Total de iterações: $i")
    return theta[end]
end

# Resolução do item a
function monte_carlo_a(M, γ, n, theta0)
    d = Gamma(theta0, 1)
    function mc()
        x = rand(d, n)
        EMV = newton_raphson(x)
        c2 = quantile(Normal(), γ + (1 - γ) / 2)
        I1 = EMV - c2 * sqrt(1 / (n * trigamma(EMV)))
        I2 = EMV + c2 * sqrt(1 / (n * trigamma(EMV)))
        return I1, I2
    end

    inferiores = []; superiores = []
    for _ in 1:M
        intervalo = mc()
        push!(inferiores, intervalo[1])
        push!(superiores, intervalo[2])
    end

    acertos = inferiores .<= theta0 .<= superiores
    return mean(acertos)

end

Rodando o código para o item a, temos

\(\text{Confiança para } 10000 \text{ simulações com alvo } 0.95, n = 100, \theta_0 = 100:\ 95.07\%\)

\(\text{Confiança para } 10000 \text{ simulações com alvo } 0.99, n = 100, \theta_0 = 100:\ 99.03\%\)

\(\text{Confiança para } 10000 \text{ simulações com alvo } 0.99, n = 5, \theta_0 = 75:\ 98.91\%\)

Para o item b,

# Resolução do item b
function monte_carlo_b(M, γ, n, theta0)
    d = Gamma(theta0, 1)
    g(a) = cdf(Gamma(a, 1), 950)
    g1(a) = pdf(Gamma(a, 1), 950)
    function mc()
        x = rand(d, n)
        EMV = newton_raphson(x)
        c2 = quantile(Normal(), γ + (1 - γ) / 2)
        I1 = g(EMV) - c2 * sqrt(g1(EMV)^2 / (n * trigamma(EMV)))
        I2 = g(EMV) + c2 * sqrt(g1(EMV)^2 / (n * trigamma(EMV)))
        return I1, I2
    end

    inferiores = []; superiores = []
    for _ in 1:M
        intervalo = mc()
        push!(inferiores, intervalo[1])
        push!(superiores, intervalo[2])
    end

    acertos = inferiores .<= g(theta0) .<= superiores
    return mean(acertos)

end

Rodando o código para o item b, temos

\(\text{Confiança para } 10000 \text{ simulações com alvo } 0.95, n = 100, \theta_0 = 1000:\ 94.88\%\)

\(\text{Confiança para } 10000 \text{ simulações com alvo } 0.99, n = 100, \theta_0 = 951:\ 98.61\%\)

\(\text{Confiança para } 10000 \text{ simulações com alvo } 0.99, n = 5, \theta_0 = 951:\ 94.67\%\)

Para o item c, precisaremos calcurar as derivadas: \[ \begin{aligned} g(\theta) &= \ln f_\theta(\theta) = - \ln \Gamma(\theta) + (\theta-1) \ln (\theta) - \theta \\ \Rightarrow g'(\theta) &= -\psi_1(\theta) + \ln(\theta) + \frac{\theta -1}{\theta} - 1 \\ &= -\psi_1(\theta) + \ln(\theta) - \frac{-1}{\theta} = g1(\theta) \end{aligned} \]

# Resolução do item c
function monte_carlo_c(M, γ, n, theta0)
    d = Gamma(theta0, 1)
    g(a) = log(pdf(Gamma(a, 1), a))
    g1(a) = -digamma(a) + log(a) - 1 / a
    function mc()
        x = rand(d, n)
        EMV = newton_raphson(x)
        c2 = quantile(Normal(), γ + (1 - γ) / 2)
        I1 = g(EMV) - c2 * sqrt(g1(EMV)^2 / (n * trigamma(EMV)))
        I2 = g(EMV) + c2 * sqrt(g1(EMV)^2 / (n * trigamma(EMV)))
        return I1, I2
    end

    inferiores = []; superiores = []
    for _ in 1:M
        intervalo = mc()
        push!(inferiores, intervalo[1])
        push!(superiores, intervalo[2])
    end
    acertos = inferiores .<= g(theta0) .<= superiores
    return mean(acertos)
end

Rodando o código para o item c, temos

\(\text{Confiança para } 10000 \text{ simulações com alvo } 0.95, n = 100, \theta_0 = 100:\ 94.77\%\)

\(\text{Confiança para } 10000 \text{ simulações com alvo } 0.99, n = 100, \theta_0 = 100:\ 99.06\%\)

\(\text{Confiança para } 10000 \text{ simulações com alvo } 0.99, n = 5, \theta_0 = 75:\ 98.97\%\)

Para o item d, note que \[ g'(\theta) = \mathrm{e}^{\ln f_\theta(x)} (-\psi_1(\theta) + \ln(\theta) - 1/\theta) \]

# Resolução do item d
function monte_carlo_d(M, γ, n, theta0)
    d = Gamma(theta0, 1)
    g(a) = exp(log(pdf(Gamma(a, 1), a))) # usado a expp(log), poderia usar direto
    g1(a) = exp(log(pdf(Gamma(a, 1), a))) * (-digamma(a) + log(a) - 1 / a)
    function mc()
        x = rand(d, n)
        EMV = newton_raphson(x)
        c2 = quantile(Normal(), γ + (1 - γ) / 2)
        I1 = g(EMV) - c2 * sqrt(g1(EMV)^2 / (n * trigamma(EMV)))
        I2 = g(EMV) + c2 * sqrt(g1(EMV)^2 / (n * trigamma(EMV)))
        return I1, I2
    end

    inferiores = []; superiores = []
    for _ in 1:M
        intervalo = mc()
        push!(inferiores, intervalo[1])
        push!(superiores, intervalo[2])
    end
    acertos = inferiores .<= g(theta0) .<= superiores
    return mean(acertos)
end

Rodando o código para o item d, temos

\(\text{Confiança para } 10000 \text{ simulações com alvo } 0.95, n = 100, \theta_0 = 100:\ 95.28\%\)

\(\text{Confiança para } 10000 \text{ simulações com alvo } 0.99, n = 100, \theta_0 = 100:\ 99.02\%\)

\(\text{Confiança para } 10000 \text{ simulações com alvo } 0.99, n = 5, \theta_0 = 75:\ 98.83999999999999\%\)

39.3.2.1 Comparação entre o exato e aproximado

using Plots, LaTeXStrings, Distributions, StatsBase, Random

function monte_carlo_realxaprox()
    theta0 = 10
    M = 100_000
    nn = 100

    γ = 0.95

    g(a) = a
    g1(a) = 1

    w = 0
    cober = zeros(nn)
    cobera = zeros(nn)
    for n in 1:nn
        w += 1
        I1 = zeros(M)
        Ia1 = zeros(M)
        I2 = zeros(M)
        Ia2 = zeros(M)

        for i in 1:M
            x = rand(Exponential(1 / theta0), n)
            c1 = quantile(Chisq(2n), (1 - γ) / 2)
            c2 = quantile(Chisq(2n), γ + (1 - γ) / 2)
            I1[i] = c1 / (2 * sum(x))
            I2[i] = c2 / (2 * sum(x))

            ca2 = quantile(Normal(), γ + (1 - γ) / 2)
            EMV = 1 / mean(x)
            Ia1[i] = g(EMV) - ca2 * sqrt(g1(EMV)^2 * EMV^2 / n)
            Ia2[i] = g(EMV) + ca2 * sqrt(g1(EMV)^2 * EMV^2 / n)
        end
        cober[w] = mean(I1 .<= theta0 .<= I2)
        cobera[w] = mean(Ia1 .<= g(theta0) .<= Ia2)
    end
    preal = scatter(
        1:10:100, cober, color = "blue",
        label = "Cobertura do IC real observado",
        title = "Comparação entre ICs: Cobertura",
        ylims = (γ - 0.1, γ + 0.1),
        xlabel = "Cobertura",
        ylabel = L"n"
    )
    scatter!(1:10:100, cobera, color = "red", label = "Cobertura IC aproximado")
    hline!([γ], label = L"\gamma")
    return preal
end

39.3.2.2 Exemplo poisson

Note que0 \[ \begin{aligned} \hat{\theta}_{\mathrm{MV}}(\boldsymbol{X}_n) &= \bar{X} \\ \sqrt{n}\frac{\hat{\theta}_{\mathrm{MV}}(\boldsymbol{X}_n) - \theta}{\sqrt{\theta}} &\stackrel{\mathcal{D}}{\rightarrow} N(0,1) \\ \sqrt{n}\frac{\hat{\theta}_{\mathrm{MV}}(\boldsymbol{X}_n) - \theta}{\sqrt{\hat{\theta}_{\mathrm{MV}}(\boldsymbol{X}_n)}} &\stackrel{\mathcal{D}}{\rightarrow} N(0,1) \\ \end{aligned} \]

function monte_carlo_poiss()
    theta0 = 10
    d = Poisson(theta0)
    g(a) = a
    g1(a) = 1
    nn = 1000
    amplitude = zeros(nn)
    cober = zeros(nn)
    M = 10_000
    for n in nn
        I1 = zeros(M)
        I2 = zeros(M)
        for i in 1:M
            x = rand(d, n)
            EMV = mean(x)
            c2 = quantile(Normal(), γ + (1 - γ) / 2)
            I1[i] = EMV - c2 * sqrt(g1(EMV)^2 * EMV / n)
            I2[i] = EMV + c2 * sqrt(g1(EMV)^2 * EMV / n)
        end
        amplitude[n] = mean(I2 .- I1)
        cober[n] = mean(I1 .<= g(theta0) .<= I2)
    end
    amp = scatter(
        10:100:1000, amplitude, title = "Amplitude dos ICs", label = "",
        color = :blue, xlabel = "n"
    )
    cob = scatter(
        10:100:1000, cober, title = "Cobertura dos ICs", label = "",
        color = :red, xlabel = "n"
    )
    plt = plot(cob, amp)
    return plt
end

39.4 Regiões de confiança

Dizemos que \(\mathrm{RC}^{(a)}(g(\theta),\gamma)\) é uma região de confança para \(g(\theta)\) com coeficiente de confiança \(\gamma\) se, e somente se, \[ P_\theta(\theta \in \mathrm{RC}(g(\theta),\gamma)) = \gamma, \forall \theta \in \Theta. \]

Se \(\lim_{n\rightarrow \infty} P_\theta(\theta \in \mathrm{RC}(g(\theta),\gamma)) = \gamma, \forall \theta \in \Theta\), dizemos que \(\mathrm{RC}^{(a)}(g(\theta),\gamma)\) é uma região de confiança assintótica.

Observação

Se \[ P_\theta(\theta \in \mathrm{RC}(g(\theta),\gamma)) \geq \gamma, \] então dizemos que \(\mathrm{RC}(g(\theta),\gamma))\) tem confiança de pelo menos \(\gamma\).

Observação

\(\mathrm{RC}(g(\theta),\gamma))\) é uma região aleatória que depende apenas da amostra aleatória e do coeficiente \(\gamma\)

39.4.1 Construção da Região de Confiança

Seja \(\hat{\theta}(\boldsymbol{X}_n)\) um estimador para “\(\theta\)” tal que

  1. \[ \sqrt{n}(\hat{\theta}(\boldsymbol{X}_n) - \theta) \stackrel{\mathcal{D}}{\rightarrow} N_p(\boldsymbol{0}, V_\theta), \forall \theta \in \Theta. \]

  2. \[ \sqrt{n}V_\theta^{-1/2}(\hat{\theta}(\boldsymbol{X}_n) - \theta) \stackrel{\mathcal{D}}{\rightarrow} N_p(\boldsymbol{0}, I), \forall \theta \in \Theta. \]

  3. \[ \sqrt{n}V_{\hat{\theta}(\boldsymbol{X}_n)}^{-1/2}(\hat{\theta}(\boldsymbol{X}_n) - \theta) \stackrel{\mathcal{D}}{\rightarrow} N_p(\boldsymbol{0}, I), \forall \theta \in \Theta. \]

  4. \[ n(\hat{\theta}(\boldsymbol{X}_n) - \theta)^T V_{\hat{\theta}(\boldsymbol{X}_n)}^{-1} (\hat{\theta}(\boldsymbol{X}_n) - \theta) \stackrel{\mathcal{D}}{\rightarrow} \chi^2_p, \forall \theta \in \Theta. \]

Podemos usar 4. para construir uma RC para “\(\theta\)” assintótica: \[ \mathrm{RC}^{(a)}(\theta, \gamma) = \left\{\theta \in \Theta : W_n(\theta) \leq q_\gamma\right\}, \] em que \(q_\gamma\) satisfaz \[ P(\chi^2_p \leq q_\gamma) = \gamma. \] Assim, por construição, \[ \lim_{n\rightarrow\infty} P_\theta(\theta \in \mathrm{RC}^{(a)}(\theta, \gamma)) = \gamma, \forall \theta \in \Theta. \]

Para o caso \(p=1\), ou seja, \(\Theta \in \mathbb{R}\), note que \[ n(\hat{\theta}(\boldsymbol{X}_n) - \theta)^T V_{\hat{\theta}(\boldsymbol{X}_n)}^{-1} (\hat{\theta}(\boldsymbol{X}_n) - \theta) = \frac{n(\hat{\theta}(\boldsymbol{X}_n) - \theta)^2}{V_{\hat{\theta}(\boldsymbol{X}_n)}} = W_n(\theta). \] Portanto, \[ \begin{aligned} \mathrm{RC}(\theta,\gamma) &= \left\{ \theta \in \Theta : \frac{n(\hat{\theta}(\boldsymbol{X}_n) - \theta)^2}{V_{\hat{\theta}(\boldsymbol{X}_n)}} \leq q_\gamma \right\} \\ \iff \mathrm{RC}(\theta,\gamma) &= \left\{ \theta \in \Theta : \lvert \hat{\theta}(\boldsymbol{X}_n) - \theta \rvert \leq \sqrt{q_\gamma \cdot \frac{V_{\hat{\theta}(\boldsymbol{X}_n)}}{n}} \right\} \\ \iff \mathrm{RC}(\theta,\gamma) &= \left\{ \theta \in \Theta : -\sqrt{q_\gamma \cdot \frac{V_{\hat{\theta}(\boldsymbol{X}_n)}}{n}} \leq \hat{\theta}(\boldsymbol{X}_n) - \theta \leq \sqrt{q_\gamma \cdot \frac{V_{\hat{\theta}(\boldsymbol{X}_n)}}{n}} \right\} \\ \iff \mathrm{RC}(\theta,\gamma) &= \left\{ \theta \in \Theta : \hat{\theta}(\boldsymbol{X}_n) -\sqrt{q_\gamma} \sqrt{\frac{V_{\hat{\theta}(\boldsymbol{X}_n)}}{n}} \leq \theta \leq \hat{\theta}(\boldsymbol{X}_n) + \sqrt{q_\gamma} \sqrt{\frac{V_{\hat{\theta}(\boldsymbol{X}_n)}}{n}} \right\} \end{aligned} \] em que \(q_\gamma\) é tal que \(P(\chi^2_1 \leq q_\gamma) = \gamma\).

No caso em que \(\Theta\) é uma semireta ou intervalo, a região de confiança coincide com o intervalo de confiança.

Observação

\[ P(\chi^2_1 \leq q_\gamma) = P(-\sqrt{q_\gamma} \leq N(0,1) \leq \sqrt{q_\gamma}) = \gamma \]

Para o caso \(p=2\), ou seja, \(\Theta = \mathbb{R}^2\), a região de confiança aproximada é \[ \mathrm{RC}(\theta, \gamma) = \left\{\theta \in \Theta : W_n(\theta) \leq q_\gamma\right\}, \]