Calcule$\lambda$ en la línea de ajuste$x^\lambda$, donde$x \in [0, 1]$

Question

Problema

Me gustaría estimación de $\lambda$ en el conjunto de la línea de $x^\lambda$, donde $x \in [0, 1]$.

Tenga en cuenta que los siguientes R código genera "cóncavo" crecimiento a medida que x aumenta de $0$ a $1$.

Lambda = 1/2.42

x = rbeta(1e4, shape1=2,shape2=2)
y = x^Lambda + rnorm(1e4, sd=.1)

plot(x,y)

Trate de

Mi enfoque es construir una función de pérdida, la L2 pérdida para encontrar el óptimo $\lambda$, R código.

SumSq = function(lam) sum((y - x^lam )^2)
optimize(SumSq, c(0,1), tol=1e-4, maximum = F )

Este código supone que voy a saber $\lambda \in (0,1)$, e optimize me da bastante resultado preciso.

Pero este enfoque no puede dar ninguna estadística asintótica resultados tales como CI, que es información útil si me gustaría considerar la incertidumbre.

Hay alguna forma estándar de hacerlo?

Answer 1

Si usted sabe cómo y = f(x) puede tratar de alinear f(x) y, a continuación, utilizar una regresión. En este caso, los registros de ambos lados, de manera que:

log(y) = lambda log(x).

Lambda = 1/2.42
x = rbeta(1e4, shape1=2,shape2=2)
y = x^Lambda + rnorm(1e4, sd=.1)

plot(x,y)
reg <- lm(log(y)~ log(x))
summary(reg)
Lambda

Tenga en cuenta que el cálculo de la lambda tiene intervalos de confianza y contiene la verdadera lambda.

Nota: Este ejemplo tiene una baja/los bien portados error que tomar el registro de error no cambia el error que mucho, log(N(0,pequeñas sigma)) ~ N(0,sigma). Si el error es suficientemente amplio/wonky suficiente, puedes probar a utilizar un modelo Lineal Generalizado (glm) o robusto modelo lineal (rlm) a cuenta de que.

Answer 2

En primer lugar, le sugiero que lea un libro excelente por Ben Bolker derecho de Modelos Ecológicos y de Datos en R. Escrito por ecologistas , creo que este es uno de los mejores libros sobre la práctica de análisis de datos, independientemente de su origen, soy ingeniero y he usado un montón. Estoy seguro de que hay otras estadísticas matemáticas libro, sin embargo este libro es, con mucho, el más práctico libro que he leído. Él también tiene un paquete de BBMLE que usted puede ser que desee comprobar. El libro, a diferencia de cualquier otro va en la estimación de máxima verosimilitud, el perfil de la probabilidad y el intervalo de confianza de la estimación y todo eso.

Volviendo a tu problema, usted necesita escribir un log-verosimilitud de la función que usted está tratando de ajustarse a los datos. Es evidente que se está suponiendo que el error se distribuye normalmente como en la simulación, por lo que el Registro de probabilidad es:

$Log\ Likelihood (Lambda,sd) = -\frac{n}{2} log(sd) - \frac{n}{2} log(2\pi)- \frac{(y-x^{Lambda})^2}{2sd} $

Usted necesita para maximizar la ecuación anterior mediante una optimización de rutina, tales como optim en R. El uso de la Hessian de la matriz a partir de la optimización para evaluar la incertidumbre de la estimación del registro de probabilidad de la función es decir, del grado de inclinación o de cómo el plano de la curvatura de la función en el punto óptimo. Si la función es caro, lo que implica menos incertidumbres en el punto óptimo tendría una estrecha banda de confianza en sus estimaciones de los parámetros, por otro lado, si su piso que habría un amplio intervalo de confianza. Hesse, matriz de Información de Fisher le ayudan a calcular el error estándar y el intervalo de confianza.

Aquí es cómo usted lo hace en R:

    set.seed(8345)

Lambda = 1/2.42

x = rbeta(1e4, shape1=2,shape2=2)
y = x^Lambda + rnorm(1e4,mean = 0, sd=.1)

plot(x,y)

## Write Log Likelihood function

log.lik <- function(theta,y,x){

  Lam <- theta[1]
  sigma2 <- theta[2]

  # sample size
  n <-  length(y)

  #error
  e<-y-(x^Lam)

  #log likelihood
  logl<- -.5*n*log(2*pi)-.5*n*log(sigma2)-((t(e)%*%e)/(2*sigma2))

  return(-logl) # R optim does minimize so to maximize  multiply by -1
}

## Estimate Paramters thru maximum likelihood

max.lik <- optim(c(1,1), fn=log.lik, method = "L-BFGS-B", lower = c(0.00001,0.00001), hessian = T,y=y,x=x)

# Lambda
Lam <- max.lik$par[1]
#0.4107119

#Fisher Information MAtrix
fisher_info<-solve(max.lik$hessian)
prop_sigma<-sqrt(diag(fisher_info))

## Estimate 95% Confidence Interval
upper<-max.lik

Answer 3

Hay alguna forma estándar de hacerlo?

Si usted piensa que la siguiente es una buena aproximación (el modelo a simular a partir de un ejemplo)

$$y_i = x_i^\lambda +\epsilon_i, \qquad \epsilon_i\sim N(0,\sigma^2)$$

i.e.,

$$\log(E(y_i))=\lambda \log x_i$$

a continuación, puede utilizar glm con family = gaussian("log") que es exactamente este modelo

lambda <- 1/2.42

set.seed(49564503)
n <- 1e4
x <-  rbeta(n, shape1 = 2,shape2 = 2)
y <-  x^lambda + rnorm(n, sd = .1)

fit <- glm(y ~ log(x) - 1, family = gaussian("log"), start =  c(0, 1))
summary(fit)
#R 
#R Call:
#R glm(formula = y ~ log(x) - 1, family = gaussian("log"), start = 1)
#R 
#R Deviance Residuals: 
#R      Min        1Q    Median        3Q       Max  
#R -0.42452  -0.06767   0.00090   0.07020   0.34418  
#R 
#R Coefficients:
#R        Estimate Std. Error t value Pr(>|t|)    
#R log(x) 0.410666   0.001807   227.3   <2e-16 ***
#R ---
#R Signif. codes:  0 ‘***' 0.001 ‘**' 0.01 ‘*' 0.05 ‘.' 0.1 ‘ ' 1
#R 
#R (Dispersion parameter for gaussian family taken to be 0.01033428)
#R 
#R     Null deviance: 1057.81  on 10000  degrees of freedom
#R Residual deviance:  103.33  on  9999  degrees of freedom
#R AIC: -17341
#R 
#R Number of Fisher Scoring iterations: 5

Observe que tanto el parámetro de dispersión y el coeficiente de estimación partido como se esperaba. Ahora, usted puede hacer que el intervalo de confianza con el error estándar por encima o por el uso de confint.