Development Tip

R의 히스토그램에 정규 곡선 오버레이

yourdevel 2020. 12. 25. 10:33
반응형

R의 히스토그램에 정규 곡선 오버레이


R의 히스토그램에 정규 곡선을 오버레이하는 방법을 온라인으로 찾았지만 히스토그램의 일반적인 "주파수"y 축을 유지하고 싶습니다. 아래의 두 코드 세그먼트를보고 두 번째 코드에서 y 축이 "밀도"로 어떻게 바뀌는 지 확인하십시오. 첫 번째 플롯 에서처럼 y 축을 "주파수"로 유지하려면 어떻게해야합니까?

보너스로 : 밀도 곡선에도 SD 영역 (최대 3 SD)을 표시하고 싶습니다. 어떻게 할 수 있습니까? 나는 시도 abline했지만 선이 그래프의 맨 위로 확장되어보기 흉해 보입니다.

g = d$mydata
hist(g)

여기에 이미지 설명 입력

g = d$mydata
m<-mean(g)
std<-sqrt(var(g))
hist(g, density=20, breaks=20, prob=TRUE, 
     xlab="x-variable", ylim=c(0, 2), 
     main="normal curve over histogram")
curve(dnorm(x, mean=m, sd=std), 
      col="darkblue", lwd=2, add=TRUE, yaxt="n")

여기에 이미지 설명 입력

위 이미지에서 y 축이 "밀도"임을 확인하십시오. 나는 그것을 "주파수"로 만들고 싶습니다.


내가 찾은 좋은 쉬운 방법은 다음과 같습니다.

h <- hist(g, breaks = 10, density = 10,
          col = "lightgray", xlab = "Accuracy", main = "Overall") 
xfit <- seq(min(g), max(g), length = 40) 
yfit <- dnorm(xfit, mean = mean(g), sd = sd(g)) 
yfit <- yfit * diff(h$mids[1:2]) * length(g) 

lines(xfit, yfit, col = "black", lwd = 2)

hist개체 에서 쉽게 계산할 수있는 올바른 승수를 찾기 만하면 됩니다.

myhist <- hist(mtcars$mpg)
multiplier <- myhist$counts / myhist$density
mydensity <- density(mtcars$mpg)
mydensity$y <- mydensity$y * multiplier[1]

plot(myhist)
lines(mydensity)

여기에 이미지 설명 입력

평균에서 멀어지는 각 표준 편차에서 정규 밀도 및 선 (평균 포함)을 사용하는보다 완전한 버전 :

myhist <- hist(mtcars$mpg)
multiplier <- myhist$counts / myhist$density
mydensity <- density(mtcars$mpg)
mydensity$y <- mydensity$y * multiplier[1]

plot(myhist)
lines(mydensity)

myx <- seq(min(mtcars$mpg), max(mtcars$mpg), length.out= 100)
mymean <- mean(mtcars$mpg)
mysd <- sd(mtcars$mpg)

normal <- dnorm(x = myx, mean = mymean, sd = mysd)
lines(myx, normal * multiplier[1], col = "blue", lwd = 2)

sd_x <- seq(mymean - 3 * mysd, mymean + 3 * mysd, by = mysd)
sd_y <- dnorm(x = sd_x, mean = mymean, sd = mysd) * multiplier[1]

segments(x0 = sd_x, y0= 0, x1 = sd_x, y1 = sd_y, col = "firebrick4", lwd = 2)

이것은 앞서 언급 한 StanLe의 anwer의 구현이며 밀도를 사용할 때 그의 대답이 곡선을 생성하지 않는 경우를 수정합니다.

이것은 기존의 숨겨진 hist.default()함수를 대체하여 normalcurve매개 변수 (기본값 :) 만 추가합니다 TRUE.

처음 세 줄은 패키지 구축 을 위해 roxygen2 를 지원하는 입니다.

#' @noRd
#' @exportMethod hist.default
#' @export
hist.default <- function(x,
                         breaks = "Sturges",
                         freq = NULL,
                         include.lowest = TRUE,
                         normalcurve = TRUE,
                         right = TRUE,
                         density = NULL,
                         angle = 45,
                         col = NULL,
                         border = NULL,
                         main = paste("Histogram of", xname),
                         ylim = NULL,
                         xlab = xname,
                         ylab = NULL,
                         axes = TRUE,
                         plot = TRUE,
                         labels = FALSE,
                         warn.unused = TRUE,
                         ...)  {

  # https://stackoverflow.com/a/20078645/4575331
  xname <- paste(deparse(substitute(x), 500), collapse = "\n")

  suppressWarnings(
    h <- graphics::hist.default(
      x = x,
      breaks = breaks,
      freq = freq,
      include.lowest = include.lowest,
      right = right,
      density = density,
      angle = angle,
      col = col,
      border = border,
      main = main,
      ylim = ylim,
      xlab = xlab,
      ylab = ylab,
      axes = axes,
      plot = plot,
      labels = labels,
      warn.unused = warn.unused,
      ...
    )
  )

  if (normalcurve == TRUE & plot == TRUE) {
    x <- x[!is.na(x)]
    xfit <- seq(min(x), max(x), length = 40)
    yfit <- dnorm(xfit, mean = mean(x), sd = sd(x))
    if (isTRUE(freq) | (is.null(freq) & is.null(density))) {
      yfit <- yfit * diff(h$mids[1:2]) * length(x)
    }
    lines(xfit, yfit, col = "black", lwd = 2)
  }

  if (plot == TRUE) {
    invisible(h)
  } else {
    h
  }
}

빠른 예 :

hist(g)

여기에 이미지 설명 입력

날짜의 경우 약간 다릅니다. 참고로 :

#' @noRd
#' @exportMethod hist.Date
#' @export
hist.Date <- function(x,
                      breaks = "months",
                      format = "%b",
                      normalcurve = TRUE,
                      xlab = xname,
                      plot = TRUE,
                      freq = NULL,
                      density = NULL,
                      start.on.monday = TRUE,
                      right = TRUE,
                      ...)  {

  # https://stackoverflow.com/a/20078645/4575331
  xname <- paste(deparse(substitute(x), 500), collapse = "\n")

  suppressWarnings(
    h <- graphics:::hist.Date(
      x = x,
      breaks = breaks,
      format = format,
      freq = freq,
      density = density,
      start.on.monday = start.on.monday,
      right = right,
      xlab = xlab,
      plot = plot,
      ...
    )
  )

  if (normalcurve == TRUE & plot == TRUE) {
    x <- x[!is.na(x)]
    xfit <- seq(min(x), max(x), length = 40)
    yfit <- dnorm(xfit, mean = mean(x), sd = sd(x))
    if (isTRUE(freq) | (is.null(freq) & is.null(density))) {
      yfit <- as.double(yfit) * diff(h$mids[1:2]) * length(x)
    }
    lines(xfit, yfit, col = "black", lwd = 2)
  }

  if (plot == TRUE) {
    invisible(h)
  } else {
    h
  }
}

참조 URL : https://stackoverflow.com/questions/20078107/overlay-normal-curve-to-histogram-in-r

반응형