>source

이것은 R에서 비선형 모델을 피팅하려는 첫 시도이므로, 제발 참아주세요.

문제

nls() 를 이해하려고 노력하고 있습니다  이 오류가 발생했습니다 :

Error in nlsModel(formula, mf, start, wts): singular gradient matrix at initial parameter estimates

가설

여기서 다른 질문에서 읽은 내용에서 다음 중 하나 일 수 있습니다.

  • 모델이 불연속 적이거나
  • 내 모델이 과도하게 결정되었거나
  • 시작 매개 변수 값의 잘못된 선택

따라서이 오류를 극복하는 방법에 대한 도움을 요청하고 있습니다. 모델을 변경하고 nls() 를 계속 사용할 수 있습니까 또는 nls.lm 를 사용해야합니까?   minpack.lm 에서  다른 곳에서 읽은 패키지?

내 접근 방식

모델에 대한 자세한 내용은 다음과 같습니다.

  • 모델은 불연속 함수, 일종의계단함수 유형 (아래 그림 참조)
  • 일반적으로 모델의단계 수는 가변적이지만 특정 피팅 이벤트에 대해 고정되어 있습니다
문제를 보여주는 MWE MWE 코드에 대한 간략한 설명
  • step_fn(x, min = 0, max = 1) : 1 를 리턴하는 함수  구간 내 ( min max ] 및 0  그렇지 않으면;이름에 대해 죄송합니다. 이제는 실제로 단계 함수가 아닙니다. interval_fn()  더 적절하다고 생각합니다.
  • staircase(x, dx, dy) : step_fn() 의 요약  기능. dx  단계의 너비로 구성된 벡터입니다 (예 : max - min ). 및 dy   y 단위로 증가  각단계마다
  • staircase_formula(n = 1L) : formula 를 생성합니다  함수 staircase() 에 의해 모델링 된 모델을 나타내는 객체  ( nls() 와 함께 사용)  기능).
  • purrr 를 사용한다는 점에 유의하십시오  그리고 glue  아래 예의 패키지.
코드
step_fn <- function(x, min = 0, max = 1) {
  y <- x
  y[x > min & x <= max] <- 1
  y[x <= min] <- 0
  y[x > max] <- 0
  return(y)
}
staircase <- function(x, dx, dy) {
  max <- cumsum(dx)
  min <- c(0, max[1:(length(dx)-1)])
  step <- cumsum(dy)
  purrr::reduce(purrr::pmap(list(min, max, step), ~ ..3 * step_fn(x, min = ..1, max = ..2)), `+`)
}

staircase_formula <- function(n = 1L) {
  i <- seq_len(n)
  dx <- sprintf("dx%d", i)
  min <-
    c('0', purrr::accumulate(dx[-n], .f = ~ paste(.x, .y, sep = " + ")))
  max <- purrr::accumulate(dx, .f = ~ paste(.x, .y, sep = " + "))
  lhs <- "y"
  rhs <-
    paste(glue::glue('dy{i} * step_fn(x, min = {min}, max = {max})'),
          collapse  = " + ")
  sc_form <- as.formula(glue::glue("{lhs} ~ {rhs}")) 
  return(sc_form)
}

x <- seq(0, 10, by = 0.01)
y <- staircase(x, c(1,2,2,5), c(2,5,2,1)) + rnorm(length(x), mean = 0, sd = 0.2)
plot(x = x, y = y)
lines(x = x, y = staircase(x, dx = c(1,2,2,5), dy = c(2,5,2,1)), col="red")


my_data <- data.frame(x = x, y = y)
my_model <- staircase_formula(4)
params <- list(dx1 = 1, dx2 = 2, dx3 = 2, dx4 = 5,
               dy1 = 2, dy2 = 5, dy3 = 2, dy4 = 1)
m <- nls(formula = my_model, start = params, data = my_data)
#> Error in nlsModel(formula, mf, start, wts): singular gradient matrix at initial parameter estimates

도움을 주시면 감사하겠습니다.

  • 답변 # 1

    길이 len 의 관측치 벡터가 주어진다고 가정합니다  귀하의 예에 표시된 것과 같이, 당신은 k 를 식별하고자합니다  점프와 k  점프 크기. (또는 어쩌면 내가 당신을 오해했지만, 당신은 정말로 당신이 달성하고자하는 것을 말하지 않았습니다.) 아래에서는 로컬 검색을 사용하여 솔루션을 스케치합니다. 나는 당신의 예제 데이터로 시작합니다 :

    x <- seq(0, 10, by = 0.01)
    y <- staircase(x,
                   c(1,2,2,5),
                   c(2,5,2,1)) + rnorm(length(x), mean = 0, sd = 0.2)
    
    

    솔루션은 점프의위치크기목록입니다. 예를 들어 20 번 점프 할 때 변수를 정의하는 것은 번거로워 지므로 벡터를 사용하여 이러한 데이터를 저장합니다.

    예제 (임의) 솔루션 :

    k <- 5   ## number of jumps
    len <- length(x)
    sol <- list(position = sample(len, size = k),
                size = runif(k))
    ## $position
    ## [1]  89 236 859 885 730
    ## 
    ## $size
    ## [1] 0.2377453 0.2108495 0.3404345 0.4626004 0.6944078
    
    

    솔루션의 품질을 계산하려면 목적 함수가 필요합니다. 또한 간단한 도우미 함수 stairs 를 정의합니다. 목적 함수에서 사용됩니다. 목적 함수 abs_diff  적합치 (솔루션에 의해 정의 된)와 y 의 평균 절대 차이를 계산합니다. .

    stairs <- function(len, position, size) {
        ans <- numeric(len)
        ans[position] <- size
        cumsum(ans)
    }
    abs_diff <- function(sol, y, stairs, ...) {
        yy <- stairs(length(y), sol$position, sol$size)
        sum(abs(y - yy))/length(y)
    }
    
    

    이제 로컬 검색의 핵심 구성 요소 인 솔루션을 발전시키는 데 사용되는 인접 기능이 제공됩니다. 이웃 함수는 해결책을 취하고 약간 변경합니다. 여기에서위치또는크기를 선택하고 약간 수정합니다

    neighbour <- function(sol, len, ...) {
        p <- sol$position
        s <- sol$size
        if (runif(1) > 0.5) {
            ## either move one of the positions ...
            i <- sample.int(length(p),  size = 1)
            p[i] <- p[i] + sample(-25:25, size = 1)
            p[i] <- min(max(1, p[i]), len)        
        } else {
            ## ... or change a jump size
            i <- sample.int(length(s), size = 1)
            s[i] <- s[i] + runif(1, min = -s[i], max = 1)
        }
        list(position = p, size = s)
    }
    
    

    전화 예 : 새 솔루션의 첫 번째 점프 크기가 ​​변경되었습니다.

    ## > sol
    ## $position
    ## [1]  89 236 859 885 730
    ## 
    ## $size
    ## [1] 0.2377453 0.2108495 0.3404345 0.4626004 0.6944078
    ## 
    ## > neighbour(sol, len)
    ## $position
    ## [1]  89 236 859 885 730
    ## 
    ## $size
    ## [1] 0.2127044 0.2108495 0.3404345 0.4626004 0.6944078
    
    

    로컬 검색을 계속 실행합니다.

    library("NMOF")
    sol.ls <- LSopt(abs_diff,
                    list(x0 = sol, nI = 50000, neighbour = neighbour),
                    stairs = stairs,
                    len = len,
                    y = y)
    
    

    해를 구할 수 있습니다. 적합 선은 파란색으로 표시됩니다.

    plot(x, y)
    lines(x, stairs(len, sol.ls$xbest$position, sol.ls$xbest$size),
          col = "blue", type = "S")
    
    

  • 답변 # 2

    대신 DE를 사용해보십시오 :

    library(NMOF)
     yf= function(params,x){
       dx1 = params[1]; dx2 = params[2]; dx3 = params[3]; dx4 = params[4];
       dy1 = params[5]; dy2 = params[6]; dy3 = params[7]; dy4 = params[8]
       dy1 * step_fn(x, min = 0, max = dx1) + dy2 * step_fn(x, min = dx1, 
                   max = dx1 + dx2) + dy3 * step_fn(x, min = dx1 + dx2, max = dx1 + 
                   dx2 + dx3) + dy4 * step_fn(x, min = dx1 + dx2 + dx3, max = dx1 + 
                   dx2 + dx3 + dx4)
     }
     algo1 <- list(printBar = FALSE,
                   nP  = 200L,
                   nG  = 1000L,
                   F   = 0.50,
                   CR  = 0.99,
                   min = c(0,1,1,4,1,4,1,0),
                   max = c(2,3,3,6,3,6,3,2))
     OF2 <- function(Param, data) { #Param=paramsj data=data2
       x <- data$x
       y <- data$y
       ye <- data$model(Param,x)
       aux <- y - ye; aux <- sum(aux^2)
       if (is.na(aux)) aux <- 1e10
       aux
     }
     data5 <- list(x = x, y = y,  model = yf, ww = 1)
     system.time(sol5 <- DEopt(OF = OF2, algo = algo1, data = data5))
     sol5$xbest
     OF2(sol5$xbest,data5)
     plot(x,y)
     lines(data5$x,data5$model(sol5$xbest, data5$x),col=7,lwd=2)
    #>  sol5$xbest
    #[1]   1.106396  12.719182  -9.574088  18.017527   3.366852   8.721374 -19.879474   1.090023
    #>  OF2(sol5$xbest,data5)
    #[1] 1000.424
    
    

  • 이전 c# - 요청시 다음 장면이로드되지 않음
  • 다음 assembly - DosBox 문자 속성을 수정하는 방법?