Resolvendo uma equação integral simples por amostragem aleatória

8

Seja uma função não negativa. Eu estou interessada em encontrar , tal que A ressalva : tudo o que posso fazer é provar em pontos em . Posso, no entanto, escolher os locais onde eu amostra aleatoriamente, se eu assim o desejarem. z [ 0 , 1 ] z 0 f ( x )fz[0,1]f [ 0 , 1 ] f

0zf(x)dx=1201f(x)dx
f[0,1]f

Questões:

  1. É possível obter uma estimativa imparcial de após muitas amostras finitas? Em caso afirmativo, qual é a menor variação possível dessa estimativa após amostras?kzk
  2. Caso contrário, quais procedimentos estão disponíveis para estimar e quais são seus tempos de convergência associados.z

Conforme apontado por Douglas Zare nos comentários, isso pode ser muito difícil se a função for próxima de zero ou muito grande. Felizmente, a função para a qual eu preciso usar isso é delimitada por cima e por baixo, então vamos supor que . Além disso, também podemos assumir que é Lipschitz ou mesmo diferenciável, se isso ajudar.f1f(x)2f

Robinson
fonte
1
Se você não tiver mais informações, poderá ter um comportamento muito ruim. Imagine que é entre e ePequenas alterações em farão a mediana saltar de abaixo de para acima de . 0 1 / 3 2 / 3 1 / 3 0 f ( x ) d x 1 / 2. f 1 / 3 2 / 3f01/32/301/3f(x) dx1/2.f1/32/3
Douglas Zare
@robinson Você poderia fornecer mais informações sobre ? Ou você está interessado em resolver o problema para qualquer densidade ? fff
@DouglasZare - Obrigado pelo comentário; veja minha edição.
Robinson
@ Procrastinator - editei a pergunta com um pouco mais de informação.
robinson
3
(+1) Para a atualização. Dividindo o lado esquerdo pela direita, pode-se ver que isso se reduz a encontrar a mediana de uma distribuição de probabilidade desconhecida suportada em . [0,1]
cardeal

Respostas:

6

Como o cardeal apontou em seu comentário, sua pergunta pode ser reafirmada da seguinte forma.

Por álgebra simples, a equação integral pode ser reescrita como na qual é a função de densidade de probabilidade definida como g g ( x ) = f ( x )

0zg(x)dx=12,
g
g(x)=f(x)01f(t)dt.

Seja uma variável aleatória com densidade . Por definição, , portanto sua equação integral é equivalente a o que significa que o seu problema pode ser declarado como:XgP{Xz}=0zg(x)dx

P{Xz}=12,

"Seja uma variável aleatória com densidade . Encontre a mediana de ".XgX

Para estimar a mediana de , use qualquer método de simulação para desenhar uma amostra dos valores de e tome como estimativa a mediana da amostra.XX

Uma possibilidade é usar o algoritmo Metropolis-Hastings para obter uma amostra de pontos com a distribuição desejada. Devido à expressão da probabilidade de aceitação no algoritmo Metropolis-Hastings, não precisamos saber o valor da constante de normalização da densidade . Portanto, não precisamos fazer essa integração.01f(t)dtg

O código abaixo usa uma forma particularmente simples do algoritmo Metropolis-Hastings conhecido como Indepence Sampler, que usa uma proposta cuja distribuição não depende do valor atual da cadeia. Eu usei propostas uniformes independentes. Para comparação, o script gera o mínimo de Monte Carlo e o resultado encontrado com a otimização padrão. Os pontos de amostra são armazenados no vetor chain, mas descartamos os primeiros pontos que formam o chamado período de "queima em" da simulação.10000

BURN_IN = 10000
DRAWS   = 100000

f = function(x) exp(sin(x))

chain = numeric(BURN_IN + DRAWS)

x = 1/2

for (i in 1:(BURN_IN + DRAWS)) {
    y = runif(1) # proposal
    if (runif(1) < min(1, f(y)/f(x))) x = y
    chain[i] = x
}

x_min = median(chain[BURN_IN : (BURN_IN + DRAWS)])

cat("Metropolis minimum found at", x_min, "\n\n")

# MONTE CARLO ENDS HERE. The integrations bellow are just to check the results.

A = integrate(f, 0, 1)$value

F = function(x) (abs(integrate(f, 0, x)$value - A/2))

cat("Optimize minimum found at", optimize(F, c(0, 1))$minimum, "\n")

Aqui estão alguns resultados:

Metropolis minimum found at 0.6005409 
Optimize minimum found at 0.601365

Esse código é apenas um ponto de partida para o que você realmente precisa. Portanto, use com cuidado.

zen
fonte
Obrigado pela sua resposta. Eu não sei R, então não tenho certeza de como analisar o que você está fazendo. Você poderia indicar em palavras / fórmulas seu procedimento? Obrigado. Em particular, estou me perguntando se você está respeitando a restrição de que a única coisa que pode fazer é avaliar f - você não tem permissão, por exemplo, para integrar (embora você possa certamente formar aproximações de Monte-Carlo para integrais com base em avaliações aleatórias). f
Robinson
Sim, estou apenas avaliando para obter a estimativa de Monte Carlo. f
Zen
O código é apenas um exemplo. A sintaxe R é semelhante a outros idiomas. Alguma declaração em particular que você não entende? Confira a página da Wikipedia sobre o algoritmo Metropolis-Hastings. Obviamente, a ideia geral é mais importante. Você pode obter amostras do usando qualquer método disponível. f/f
Zen
Você fez um curso introdutório sobre processos estocásticos, cobrindo cadeias de Markov com tempo discreto?
Zen
1
BTW: Procrastinadores do mundo, uni-vos! Mas hoje não ...
Zen
3

A qualidade da aproximação integral, pelo menos no caso tão simples quanto 1D, é dada por (Teorema 2.10 em Niederreiter (1992) ): onde é o módulo de continuidade da função (relacionado à variação total e facilmente expressável para as funções de Lipshitz) e é a discrepância (extrema) ou a diferença máxima entre a fração de pela sequência

|1Nn=1Nf(xn)01f(u)du|ω(f;DN(x1,,xN))
ω(f;t)=sup{|f(u)f(v)|:u,v[0,1],|uv|t,t>0}
DN(x1,,xN)=supu|1Nn1{xn[0,u)}u|=12N+maxn|xn2n12N|
x1,,xNde um intervalo semiaberto e sua medida de Lebesgue . A primeira expressão é a definição e a segunda expressão é a propriedade das seqüências 1D em (Teorema 2.6 no mesmo livro).[0,u)u[0,1]

Então, obviamente, para minimizar o erro na aproximação integral, pelo menos no RHS da sua equação, você precisa tomar . Dane-se as avaliações aleatórias, elas correm o risco de ter uma lacuna aleatória em uma característica importante da função.xn=(2n1)/2N

Uma grande desvantagem dessa abordagem é que você precisa se comprometer com um valor de para produzir essa sequência uniformemente distribuída. Se você não estiver satisfeito com a qualidade da aproximação que ela fornece, tudo o que você pode fazer é dobrar o valor de e atingir todos os pontos médios dos intervalos criados anteriormente.NN

Se você deseja ter uma solução em que possa aumentar o número de pontos mais gradualmente, continue lendo esse livro e aprenda sobre as seqüências de van der Corput e inversões radicais. Veja Seqüências de baixa discrepância na Wikipedia, ele fornece todos os detalhes.

Atualização: para resolver para , defina a soma parcial Encontre modo que e interpole para encontrar Essa interpolação assume que é contínuo. Se adicionalmente for duas vezes diferenciável, essa aproximação integrará a expansão de segunda ordem para incorporar e e resolver uma equação cúbica para .z

Sk=1Nn=1kf(2n12N).
kzN=2k-1
Sk12SN<Sk+1,
f()f()Sk-1Sk+2z
zN=2k12N+SN/2SkN(Sk+1Sk).
f()f()Sk1Sk+2z
StasK
fonte
1
Eu gosto da essência disso. Eu acho que seria útil tornar mais explícita a estratégia que você propõe para resolver a questão do OP. Atualmente, a resposta diz (para mim) principalmente como se estivesse abordando como calcular o RHS da equação na pergunta.
cardeal
1
(+1) Boa atualização. pode ser simplesmente visto como uma aproximação da soma de Riemann à integral, onde usamos o valor de no ponto médio de cada intervalo definido pela partição, em vez do ponto final esquerdo ou direito. :-) fSNf
cardeal
Sim; é interessante, porém, que essa soma de Riemann tenha essa justificativa de otimização.
Stask