Ajustando uma distribuição para corresponder aos pontos conhecidos no CDF

Recentemente, encontrei uma situação em que conheço alguns pontos de probabilidade na cauda de uma distribuição e quero "ajustar" uma distribuição que passa por esses pontos na cauda. Sei que isso é confuso e não muito preciso e está atormentado por questões conceituais. No entanto, confie em mim que realmente quero fazer isso.

Tão efetivamente conheço alguns pontos na cauda do CDF por xserem os valores e ya probabilidade desse valor ou menor. Aqui está o código R para ilustrar meus dados:

x <- c(0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85)
y <- c(0.0666666666666667, 0.0625, 0.0659340659340659, 0.0563106796116505, 
       0.0305676855895196, 0.0436953807740325, 0.0267459138187221)

Em seguida, crio uma função para minimizar o erro entre meus dados e um CDF de distribuição beta usando pbeta. Eu uso o SSE como uma métrica adequada e minimizo isso com -sum. Eu dou um palpite inicial como o primeiro parâmetro optimde, (9, .8)embora eu tenha tentado isso com suposições diferentes e sempre obtenho o mesmo resultado. O ponto de partida que eu acho que vem vem da elaboração manual de parâmetros manuais que parecem próximos.

# function to optomize with optim
beta_func <- function(par, x) -sum( (pbeta( x, par[1], par[2]) - y)**2 ) 
out <- optim(c(9,.8), beta_func, lower=c(1,.5), upper=c(200,200), method="L-BFGS-B", x=x)

out <- out$par
print(out)
#> [1]  0.90000 23.40294

Abaixo, eu gráfico a distribuição beta 'otimizada' em vermelho, meus dados reais em azul e uma mão aprimorou o palpite inicial dos parâmetros beta em preto.

plot(function(x) pbeta(x, shape1=out[1], shape2=out[2] ), 0, 1.5, col='red')
plot(function(x) pbeta(x, 9,.8), 0, 1.5, col='black', add=TRUE)
lines(x,y, col='blue')

Não posso grunhir o que está acontecendo optimpara dar uma solução pior que o meu palpite inicial. Calculei o SSE para meu palpite inicial versus a optimsolução e parece que meu palpite tem um -SSE muito maior:

# my guess
-sum( (pbeta( x, 9, .8) - y)**2)
#> [1] -0.03493344

# optim's output
-sum( (pbeta( x, .9, 23) - y)**2)
#> [1] -6.314587

Usando a história do passado como meu anterior bayesiano, meu palpite é que estou entendendo errado optimou alimentando-o com entradas impróprias. No entanto, eu não consigo entender o que está acontecendo. Quaisquer dicas seriam extremamente apreciadas.

Eu tentei usar o CGmétodo de otimização, mas os resultados não são significativamente diferentes e ainda não parecem tão bons quanto o meu palpite inicial.

out <- optim(c(9,.8), beta_func, method="CG", x=x)
out <- out$par
print(out)
#> [1]  2.287611 11.124736

Eu acho que você está acidentalmente tentando maximizar os erros ao quadrado. O padrão para optim () é minimizar a função, para que o sinal negativo em seu beta_func () resulte na pesquisa de um lance máx.

Se você modificar sua função assim, obterá valores mais próximos do seu palpite:

beta_func2 <- function(par, x) sum( (pbeta( x, par[1], par[2]) - y)**2 ) 
out2 <- optim(c(9,.8), beta_func2, lower=c(1,.5), upper=c(200,200), method="L-BFGS-B", x=x)
out2 <- out2$par
print(out2)
[1] 11.04296  0.50000

Você pode verificar a nova função com base nas suas observações (onde out, x e y são definidos como no seu exemplo):

plot(x,(pbeta(x,out[1],out[2])), ylim=c(-.1,1), col="red", type="l")
points(x, (pbeta(x,9,0.8)), col="black", type="l")
points(x,(pbeta(x,out2[1],out2[2])), col="orange", type="l")
lines(x,y, col='blue')
title(main="Blue observed CDF, black guesstimate, gold optimized")